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Natural  Convection  Heat  TVansfer  in 
Water  Near  its  Density  Maximum 

YIN-CHAO  YEN 


INTRODUCTION 

The  most  comprehensive  treatment  of  cold  regions  heat  transfer  processes  is  given  by  Lunardini 
(1981,  1988).  He  discusses  the  mechanisms  of  frost  heave  and  their  relation  to  soil  composition  and 
properties,  methods  of  calculating  or  approximating  the  extent  of  ground  freezing  and  thawing  with 
uniform  or  variable  thermal  properties,  and  transfer  processes  associated  with  structures  in  cold 
climates.  In  contrast,  this  monograph  deals  exclusively  with  the  phenomena  arising  from  the  density 
anomaly  of  water  that  is  formed  by  melting  ice  in  either  an  ice-water  system  or  in  single-phase  water 
that  encompasses  the  density  extremum.  The  effect  of  this  density  anomaly  on  the  onset  of  convection 
(i.e.,  thermal  instability),  the  transient  and  steady  temperature  structure  within  the  water  layer,  and  the 
associated  heat  transfer  rate  will  be  examined  and  summarized  in  various  geometrical  arrangements. 


GENERAL  DESCRIPTION 
OF  RAYLEIGH  INSTABILITY 

Forfluids  not  undergoing  phase  change  or  density  inversion,  the  magnitude  of  the  Rayleigh  number 
is  the  criterion  commonly  used  to  determine  if  natural  convection  is  significant: 


Ra  - 


(1) 


where  P  =  the  coefficient  of  expansion 
v  =  the  kinematic  viscosity 
a  =  the  thermal  diffusivity  of  the  fluid 
d  =  the  depth  of  the  fluid  layer 
g  =  the  gravitational  acceleration 
AT  =  the  temperature  difference  across  the  liquid  layer. 


P,  v  and  a  are  usually  evaluated  at  the  arithmetic  average  of  lower  and  upper  boundary  temperatures. 
For  normal  fluids  with  density  decreasing  with  increasing  temperature,  it  has  been  determined 
analytically  and  confirmed  experimentally  that  the  liquid  will  go  into  motion  at  a  Rayleigh  number  of 
about  1 708  for  both  lower  and  upper  rigid  boundaries  and  constant  temperatures. 


The  ice-water  system,  in  which  the  water  accumulates  continuously  by  melting  from  below,  forms 
in  the  lower  part  of  the  system,  and  is  subject  to  a  negative  temperature  gradient,  resembles  the  classical 
Rayleigh  problem.  However,  the  melting  system  differs  fundamentally  in  a  number  of  ways  from  that 
of  the  classical  Rayleigh  problem:  1)  for  an  ice-water  system,  water  may  not  act  as  a  normal  fluid 
because  of  its  peculiar  properties  of  density  inversion  at  4°C;  therefore,  ideally,  a  stable  and  unstable 
region  may  exist  simultaneously  in  the  water  layer,  2)  the  upper  boundary  (water/ice  interface)  is 
moving  upward  as  melting  progresses  so  that  it  can  be  considered  neither  a  rigid  nor  a  free  boundary, 
and  3)  though  both  boundaries  are  maintained  at  fixed  temperatures,  the  upper  boundary  is  moving; 
therefore,  a  time-dependent  temperature  gradient  exists  across  the  layer. 

For  the  case  of  the  ice-water  system,  with  the  water  layer  formed  by  melting  ice  from  above,  the 
water  layer  will  be  situated  over  the  remaining  ice.  For  normal  fluids,  the  system  will  always  be  stable. 
However,  in  the  case  of  water  near  the  water/ice  interface,  there  exists  a  region  with  a  positive  density 
and  temperature  gradient  directed  away  from  the  interface.  This  is  a  system  equivalent  to  the  case  of 
a  water  layer  heated  above  (£  4°C)  and  cooled  below  at  0°C  in  which  a  stable  region  will  sit  on  an 
unstable  region.  A  convective  motion  is  thus  created  in  the  layer  for  both  cases,  involving  either  phase 
change  (i.e.,  in  a  ice-water  system)  or  a  single  phase  (i.e.,  in  water  alone),  either  when  1)  the  lower 
is  maintained  at  8°C  (i.e.,  heating  from  below)  and  the  upper  boundary  is  at  0°C  (equivalent  to  the  case 
of  melting  from  below),  or  2)  the  upper  boundary  is  maintained  at  8°C  and  the  lower  is  at  0°C 
(equivalent  to  melting  from  above).  When  the  value  of  P  in  eq  1  is  evaluated  at  8/2 = 4°C  [P  =  1/p  (3p / 
9T)4oc =0],  one  may  face  the  contradictoiy  situation  of  observing  convection  with /fa =0,  which  violates 
the  classical  stability  criterion. 


EFFECT  OF  4°C  WATER 
ON  ONSET  OF  CONVECTION 

In  a  confined  horizontal  layer 

When  a  layer  of  liquid  whose  density  decreases  monotonically  with  the  increase  of  temperature  is 
subject  to  an  adverse  temperature  gradient  (say  Tx  >  T2  where  subscripts  1  and  2  refer  to  lower  and 
upper  surfaces,  respectively)  the  system  is  potentially  unstable  because  it  is  top-heavy.  The  onset  of 
convection  is  indicated  if  the  Rayleigh  number  exceeds  its  critical  value.  For  normal  fluid,  the  Rayleigh 
number  is  defined  as  in  eq  1 . 

A  somewhat  complicated  situation  arises  if  the  liquid  possesses  a  maximum  density  within  the 
temperature  range  between  Tx  and  T2  (such  a  situation  could  exist  either  in  a  water  layer  formed  by 
melting  ice,  i.e.,  ice-water  system,  or  in  a  water  layer  by  maintaining  the  boundary  layer  temperatures 
in  such  a  manner  that  a  maximum  density  will  exist  within  the  layer).  In  such  a  case,  the  liquid  density 
increases  upward  from  the  lower  surface  until  it  reaches  a  maximum  and  decreases  afterward. 
Therefore,  only  part  of  the  liquid  layer  is  potentially  unstable.  Furthermore,  the  onset  of  convection 
is  possible  with  both  heating  from  below  and  above. 

For  liquids  with  a  density-temperature  relationship  expressed  by 

P=  ftn  [l-Y  (r-Tn,)2]  (2) 

Veronis  (1963),  Debler  (1966)  and  Tien  (1968)  all  have  found  that  the  critical  Rayleigh  number  is 
dependent  upon  the  temperature  difference  ratio  defined  as  A  =  ( Tx-Tm)/(Tx-T2 )  with  the  Rayleigh 
number  given  by 

Ra  =  2Tfri-rm)g  M)j  (3) 
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where  y  is  the  temperature  coefficient  of  eq  2,  and  Tm  is  the  temperature  corresponding  to  maximum 
density. 

The  major  limitation  of  these  earlier  works  is  that  the  particular  temperature-density  relationship 
used  has  only  a  limited  range  of  applicability.  For  the  case  of  water,  the  representation  of  the 
temperature-density  relation  by  a  parabolic  expression  is  only  valid  for  the  temperature  range  from 
0°  to  8°C.  The  other  limitation  was  that  these  investigations  were  restricted  to  the  special  case  of  two 
rigid  boundary  surfaces. 

Sun  et  al.  ( 1 969)  extended  these  earlier  works  by  representing  the  density-temperature  relationship 
with 


P  =  Pm  [l-Yi  ( T-Tm )2  -  y2  (r-  Tm)3]  (4) 

which  was  found  to  be  valid  for  water  temperature  from  0°  to  30°C.  They  also  broadened  their  in¬ 
vestigation  by  including  the  rigid-free  hydrodynamic  boundary,  in  addition  to  rigid-rigid  boundary. 
Using  the  classical  stability  analysis,  they  derived  a  new  Rayleigh  number  as 


Ra  = 


2yiAg{AT)2  d2  (l  +  |  12.  A  AT) 
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in  terms  of  thermal  parameters 
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The  critical  Rayleigh  number  has  been  computed  for  the  range  of  -4.25  <  X,  <  -0.5  and  -1 .4  <  A^ 
<  1 .6.  For  the  special  case  A^ = 0  (equivalent  only  to  the  parabolic  density-temperature  relation,  where 
y,  is  only  a  function  of  the  temperature  difference  ratio  A  since  y2  =  0),  the  critical  Rayleigh  number 
for  both  rigid-free  and  rigid-rigid  boundary  conditions  is  shown  in  Figure  1 .  The  effect  of  A^  on  Rac 
is  shown  in  Figures  2  and  3,  respectively,  for  rigid-free  and  rigid-rigid  boundaries. 

It  can  be  seen  from  Figure  1,  that  forXj  >-l  .75,  the  values  of  (Rac)\^.0 are  nearly  equal  to  each  other 
for  the  conditions  of  rigid-free  and  rigid-rigid  boundaries. 

Chandrasekhar  (1961)  derived,  for  the  use  of  rigid-rigid  boundary,  the  asymptotic  expression  for 
Rac  as 


Rac  ~  1,186  ji,)4  (8) 

in  terms  of  the  work  reported  by  Sun  et  al.  (1969).  Equation  8  then  becomes 

Rac  ~  1,186  (A,])4,  for  A.2  =  0 .  (9) 
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Though  the  above  approximation  is  obtained  for  the  rigid-rigid  surfaces,  this  is  expected  to  hold  true 
for  rigid-free  surfaces  and  for  the  case  of  melting  from  below.  There  is  no  proof  in  the  case  of  melting 
from  above,  where  the  unstable  region  is  above  the  moving  water/ice  interface.  Fora  layer  of  maximum 
density  fluid  subject  to  temperatures  from  above  and  below,  which  are  on  either  side  of  the  maximum 
density  temperature,  the  density  of  the  fluid  first  increases  upward  from  the  lower  surface  to  a 
maximum  then  decreases.  The  fluid  layer,  therefore,  consists  of  a  potentially  unstable  layer  (with  a 
height  of  Ad)  and  a  stable  layer.  The  upper  boundary  effect  (i.e.,  one  with  free  or  rigid  upper  surfaces) 
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Figure  3.  (Ra^/fRa^.^  vs  ^  with  X;  as  pa-  Figure  4.  Rac  as  function  ofTj  or  T2  (after  Sun  et  al. 

rameter  for  rigid-rigid  case  (after  Sun  et  al.  1969). 

1969). 

is  felt  only  indirectly  by  the  potentially  unstable  layer.  This  boundary  effect  would  become  less 
significant  as  the  relative  thickness  of  the  stable  layertothat  of  the  potentially  unstable  layer  increases, 
which  corresponds  to  a  decreasing  value  of  A  or  an  increasing  value  of -X,  (see  Fig.  1).  For  the  limiting 
case  where  A  -»  0,  Tt  =  Tm,  and  the  fluid  layer  is  always  stable.  Figure  4  completely  supports  this 
argument.  Although  the  curve  for  the  rigid-rigid  case  differs  from  that  of  rigid-free  case  for  small 
values  of  (-A.,),  both  curves  approach  the  same  asymptotic  condition  for  (-X,)  >  3.0  and  become 
indistinguishable.  The  asymptotic  expression  is  found  to  be 

(*«c)k-o  ~  U77  (X,)4.  (10) 

This  is  essentially  the  same  as  eq  9,  with  a  less  than  1%  difference  in  coefficient  which,  in  all 
probability,  can  be  attributed  to  errors  introduced  in  the  numerical  computations. 

Experimental  verification 

Yen  ( 1 968 ),  Yen  and  Galea  ( 1 969),  Sun  et  al .  ( 1 969)  and  Yen  ( 1 980)  verified  the  above  findings  by 
performing  melting  experiments  in  both  melting  modes  ( i.e.,  melting  from  below  and  above).  In  both 
cases,  bubble-free  ice  was  fabricated,  and  the  melting  system  was  maintained  airtight  (any  trace  of  air 
trapped  in  the  water/ice  interface  would  drastically  affect  the  melting  mechanism).  In  these  experiments 
the  onset  of  convection  was  found  either  by  determining  the  inflection  point  on  the  melting  front  vs 
time  plot  or  by  determining  the  time  where  the  steep  jump  in  the  temperature  gradient  of  the  stable 
region  occurred.  If  the  values  of  T,  (T,  >  0°C,  melting  from  below)  and  T2  (J2  >  0°C  melting  from 
above)  and  the  critical  layer  depths  determined  by  either  of  the  mentioned  methods  (i.e.,  inflection 
point  or  sudden  jump  in  the  temperature  gradient)  are  known,  the  experimental  critical  Rayleigh 
number  can  be  found  from 
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Figure5.  Principal  stability  diagram(after\ierker 
et  al.  1979). 
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The  experimental  values  of  Race  from  eq  1 1  a  and  eq  1 1  b  are  compared  with  Raa  values  obtained  by 
using  A.,  and  Rvalues  from  eq  12and  13  inconjunction  with  Figures  1  and  3.  Hie  comparison  is  shown 
in  Figure  4,  Within  the  experimental  error  in  determining  the  critical  layer  depth  dc  (which  is  exem¬ 
plified  by  its  presence  as  d2  in  eq  1  la  and  1  lb),  it  can  be  concluded  that  Race  and  Racl  are  nearly  in 
complete  agreement. 

The  other  comprehensive  analysis  of  the  onset  of  convection  (i.e.,  critical  Rayleigh  number)  in  a 
water  layer  encompassing  the  density  anomaly  was  reported  by  Mericer  et  al.  (1979).  In  their  study, 
the  density-temperature  relation  of  water  was  approximated  by  three  different  polynomials  having 
two,  three  and  five  terms.  Linear  stability  analysis  was  used  and  the  resulting  perturbation  equations 
were  solved  with  the  aid  of  Galerkin’s  method.  Figure  5  is  a  general  diagram  of  regions,  in  which  the 
fluid  is  purely  stable,  unstable  or  only  unstable  at  the  lower  layer,  according  to  the  imposed  boundary 
temperatures  of  T ,  and  T2.  P ,  and  P2  are  the  coefficients  of  tnermal  expansion  evaluated  at  T,  and  T2 
respectively,  and  the  values  of  Rat  and  Ra2  are  the  corresponding  Rayleigh  numbers  defined  as  in  eq 
1  with  AT  =  T|  -  T2.  Ral  is  always  defined  as  positive  if  the  layer  is  unstably  stratified,  whereas  Ra2 
changes  sign.  Therefore  Rai  was  used  by  Merker  et  al.  (1979)  to  describe  the  layer  stability. 

Figure  6  shows  the  critical  Rayleigh  number  variation  with  7",  in  the  case  of  constant  boundary 
temperature  (i.e.,  Tw  =  constant)  with  the  density-temperature  relation  represented  by  a  polynomial 
of  order  p = 5  and  n = 7  (number  of  Galerkin  terms).  Figure  7  shows  the  case  of  constant  wall  heat  flux 
density  (i  .e.,  qw = constant).  Both  figures  show  that  in  the  region  below  the  isotherm  7*2 = 4°C,  the  water 
layer  has  a  density  profile  with  no  maximum  value;  i.e.,  the  layer  is  completely  unstably  stratified.  For 
this  case,  the  bending  of  the  density  profile  is  weak,  and  accordingly,  the  effect  on  the  critical  Rayleigh 
number  is  small.  The  calculated  Rac  values  are  between  1708  and  3600  for  Tw  =  constant  and  below 
720  and  1 600  for  qw = constant.  The  region  above  the  isotherm  T2 = 4°C  refers  to  adensity  profile  with 
maximum  value;  i.e.,  only  a  part  of  the  layer  is  unstably  stratified.  The  nonlinearity  of  the  density 
profile  is  strong  and  the  effect  on  Rac  is  great.  The  Rac  values  are  significantly  larger  than  those 
obtained  from  the  classical  problem.  The  usage  of  Figure  6  or  7  can  be  delineated  by  considering  the 
following  cases: 

1 .  If  the  temperature  T2  <  4°C  (say  T2  =  0°C)  then  cooling  this  water  layer  with  Tt  <  0°C  from 
below  results  in  an  unstable  stratification  where  the  bending  of  the  profile  is  weak.  The  water 
layer  is  stable  if  heated  from  below  with  temperature  0  £  Tj  £  4°C  and  it  becomes  partially 
unstable  if  T,  exceeds  4°C.  The  density  profile  then  includes  the  maximum  density  and  the 
bending  of  the  profile  is  strong. 


Figure  6.  Critical  Rayleigh  number  Rac  as 
function  of  T ,  with  T2  as  a  parameter for  7W 
( wall  temperature) = constant  (after  Merker 
etal.  1979). 
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Figure  7.  Critic  al  Rayleigh  number  Rac.t«  func¬ 
tion  of  T  j  withT^asa  parameter  for  qH,(  constant 
wall  flux)  =  constant  ( after  Merker  et  al.  1979). 


2.  If  T2  is  kept  at  7"2>  4°C  (say  72=8°C)  then  cooling  this  water  from  below  with  T ( <  4°C  results 
in  the  upper  part  between  the  upper  boundary  and  the  4°C  isotherm  remaining  stable,  whereas 
the  lower  part  becomes  unstable.  The  strong  bending  of  the  profile  affects  the  Rayleigh  number 
considerably.  The  layer  is  stable  if  cooled  with  temperature  T,  between  4  and  8°C,  and  is  unstably 
stratified  if  heated  from  below  with  Ti  >  8°C. 

For  constant  wall  temperature  and  rigid  boundary  conditions,  this  study  should  be  comparable  to  that 
of  Sun  et  al.  (1969)  even  though  the  Rayleigh  number  is  defined  here  classically  (except  that  p  is 
evaluated  at  the  lower  boundary  temperature  or  T ,)  and  that  of  Sun  et  al.  (1969)  is  defined  by  eq  5. 
ForT  |  =  10°C  and  T2 = 0°C  (i.e.,  the  case  of  melting  from  below)  the  values  of  A.,  and  A^  are  1 .527  and 
0.2 1 7,  respectively,  based  on  eq  1 2b  and  1 3b.  In  Figures  1  and  3,  the  Rac  value  was  found  to  be  -8,200. 
Under  the  same  conditions  of  T,  and  T2,  based  on  Figure  6,  Rac  is  found  to  be  -8, 1 50.  For  the  case  of 
melting  from  above  (i.e.,  for  the  case  7’1  =  0°C),  and  T2=  10.5°C,  based  on  eq  12a  and  13a,  the  values 
of  X.j  =  1 .644  and  A^  =  0.300  were  obtained  with  the  aid  of  Figures  1  and  3,  Rac  is  found  to  be  Rac  = 
48,200,  which  compares  remarkably  well  with  Rac  =  47,520  estimated  from  Figure  6.  It  is  clear  that 
regardless  of  how  you  define  the  Rayleigh  number  (as  long  as  done  appropriately),  the  equations 
will  provide  the  same  Rac  for  the  same  thermal  boundary  conditions. 

Legros  et  al.  ( 1 974)  were  the  first  ones  to  investigate  systematically  the  effect  of  the  ratio  d/h  (where 
d  is  total  layer  depth,  h  is  height  of  the  unstable  layer)  on  the  onset  of  convection  of  a  water  layer  with 


Table  1.  Summary  of  experimental  parameters  and  calculated  results 
(after  Legros  et  al.  1974). 
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AT 
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0.54 

4.67 

3.17 

1.50 

1.47 

10.80 

4.32 

4.67 

-0.23 

6.62 

0.86 

4.67 

3.27 

1.40 

1.42 

10.46 

5.18 

4.67 

-0.96 

5.75 

1.31 

4.67 

3.39 

1.28 

1.38 

9.84 

6.51 

4.67 

-2.47 

4.75 

1.67 

4.67 

3.51 

1.16 

1.33 

9.40 

9.05 

4.67 

-6.94 

3.40 

2.17 

4.67 

3.70 

0.97 

1.27 

8.80 

13.33 

4.67 

-19.90 

2.19 

2.56 

4.67 

3.85 

0.82 

1.21 

8.24 

16.84 

4.67 

-34.66 

1.73 

3.40 

4.67 

4.29 

0.38 

1.09 

7.38 

20.16 

4.67 

-55.90 

1.35 
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a  maximum  density  as  a  function  of  the  upper  boundary  temperature  (i.e.,  T2).  In  their  experiment,  T2 
was  varied  from  0.54°  to  20. 1 6°C,  and  the  layer  depth  was  fixed  at  0.467  cm.  The  critical  temperature 
difference  (ATce)  was  determined  by  observing  the  temperature  difference  as  a  function  of  heat  flux. 
Table  1  is  a  summary  of  their  results. 

In  Table  1 ,  d-h  =  d/ATce  (4  -  T2),  (i.e.,  the  stable  layer  depth  just  prior  to  onset  of  convection).  For 
T2>  4°C,  the  measured  A Tce  values  were  compared  with  the  analytical  A  Ta  values  calculated  from  the 
classical  critical  Rayleigh  number,  i.e., 

gMZctjj  =  1708 
va 

or 

4r„  =  <1708>  <vg> . 

sP**3 

The  p  value  was  evaluated  from  a  sixth-order  polynomial:  i.e., 

p  =  aT6  +  bT5  +  CT4  +  dTi  +  eT2  +  fT  +  g  , 

where  a  =  8.2270200x1 0"13, 
b  =  6.8247601x10-' *, 
c  =  2.9114740xl0"9, 
d  =  1.2488637xl0-7, 
e  =  -  9.2582 184XKH, 

/  =  6.8355321X10"5, 
g  =  9.9984055x10"', 

and  was  found  to  be  in  excellent  agreement.  The  difference  AT^  -  A Ta  becomes  greater  when  T2  -» 
4°C,  as  expected.  For  the  experiments  of  T2  <  4°C,  the  results  were  compared  using  an  analogy 
pointed  out  by  Veronis(  1963)  and  Debler(1966)  between  theeigenvaluesfor  the  Benard  problem  with 
a  maximum  density  and  for  the  stability  of  viscous  flow  between  two  cylinders  rotating  at  almost  the 
same  angular  velocity.  This  is  shown  in  Figure  8  along  with  the  work  reported  by  Yen  (1980)  from 
melting  ice  studies.  Yen’s  data  not  only  agree  extremely  well  with  those  of  Veronis  (1963)  and  Legros 
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Figure  8.  Comparison  of  Ta and  Ra^ 
2K4 rf  vs  \  (after  Yen  1980). 
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et  al.  (1974)  but  the  data  extended  Legros  et  al.’s  range  of  k  from  about  1 .5  to  3.5.  In  this  figure,  the 
value  of  Rac  is  evaluated  from 

,  ajATfhl 
va 

where  y  is  obtained  from  eq  2  and  has  a  value  of  7.68x1  O^C-2  with  pm  =  0.999973  g/cm3,  and  v  and 
a  are  evaluated  at  ( T y+T^H.  In  Yen’s  melting  ice  work,  AT = 4°C  and  hc = 4 /T2(dc)  for  melting  above 
and  AT  =  T,  -  4°C  and  hc  =  [(T ,  —  4)/T ,  ]  dc  for  melting  from  below,  where  dc  is  the  total  layer  depth 
at  the  onset  of  convection. 

In  a  melt  layer  with  a  free  surface 

Seki  etal.  (1977)  performed  an  experimental  and  analytical  study  on  the  thermal  instability  program 
arising  from  a  horizontal  layer  of  ice  heated  from  above  under  constant  radiant  heat  flux.  Two  cases 
were  considered,  as  indicated  in  Figure  9.  In  the  case  of  (a),  the  fluid  density  in  the  layer  increases  at 
first  downward  from  the  free  surface  and  then  decreases  (equivalent  to  melting  from  above  with 
constant  upper  boundary  T2  =  constant).  Therefore,  the  fluid  layer  consists  of  potentially  stable  and 
unstable  layers.  On  the  other  hand,  in  the  case  of  (b),  the  entire  layer  is  unstable.  In  the  case  of  (a),  the 
hydrodynamic  boundary  conditions,  including  surface  tension  and  the  thermal  boundary  conditions, 
are  only  felt  indirectly  by  the  potentially  unstable  layer.  If  the  thickness  of  the  stable  layer  decreases 
relative  to  that  of  the  unstable  one,  it  is  obvious  that  these  effects  would  become  more  significant.  In 
case  of  (b),  these  effects  are  felt  directly  by  the  potentially  unstable  layer.  The  problem,  which 
incorporates  a  maximum  density  and  an  upper  free  surface,  differs  distinctly  from  other  instability 
problems. 

In  Seki  et  al.’s  (1977)  analysis,  linear  perturbation  techniques  are  used  to  derive  a  a  sixth-order 
differential  equation  and  the  series-solution  method  is  utilized  to  obtain  an  eigenvalue  equation  for 
the  case  where  the  lower  surface  is  kept  at  0°C  and  the  upper  free  surface  is  subjected  to  temperatures 
ranging  from  1 .5°  to  12°C.  Figure  10  shows  the  variation  of  the  critical  Rayleigh  number  Rac  with  T2 
for  Ma  =  0  (Afa  is  the  Marangoni  number).  Two  solid  lines  were  shown.  One  for  Bi = 0  (Biot  number 
=  hH/K)  corresponds  to  a  fixed  heat  flux  and  the  other  for  Bi  =  «>  to  a  fixed  temperature  at  the  free 
surface.  Also  shown  are  the  results  of  Sun  et  al.  (1969).  It  appears  that  Rac  values  increase  with  T2. 
However,  the  two  curves  intersect  each  other  at  ~6.2°C,  and  the  reason  for  this  occurrence  is  still  not 
understood.  (This  phenomenon  was  also  observed  in  Sun  et  al.’s  study;  see  Fig.  1 .)  The  effect  of  Bi 
on  Rac  is  demonstrated  through  the  plot  of  RaJ(Rac)Bi=Q  vs  Bi  for  various  T2  for  Afa=0  (see  Fig.  1 1 ). 
It  is  clear  that  the  effect  is  greatest  for  fixed  surface  temperature  and  smallest  for  fixed  heat  flux.  But 
the  sensitivity  to  Bi  varies  considerably;  it  decreases  with  increasing  T2,  indicating  that  the  effect 
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Figure  9.  Illustration  of  the  effect  of  maximum  density  in  a  melted 
water  layer  (after  Seki  et  al.  1977). 


10 


Figure  10.  Critical  Rayleigh  number  as  function  o/T2  (Seki  et 
al.  1977). 


becomes  less  significant  as  the  thickness  of  the  stable  layer  increases.  The  effect  of  Ma  on  the  ratio  of 
Ra(J(Rac)Mi  _o  vs  Ma  for  various  Bi  values  and  T2  =  4°C  is  shown  in  Figure  1 2.  Clearly  the  effect  of 
thermal  boundary  condition  at  the  free  surface  on  the  onset  of  free  convection  for  a  given  Ma  is  great¬ 
est  for  Bi  =  0  and  smallest  for  Bi  =  OO, 

The  dependence  of  /toc/(/?<Jc)Ma=0on  Ma  is  shown  in  Figure  13  with  Bi  =  0  and  various  T2  values. 
/?ac/(/?ac)Ma=0  increases  monotonically  with  increasing  Ma  and  is  most  sensitive  to  Ma  in  the  range 
of  100-1000.  The  value  of  T2markedly  influences  the  effect  of  Ma  on  Rac.  For  T2>  4°C,  the  thickness 
of  the  potentially  stable  layer  in  the  melted  water  layer  increases  as  T2  increases;  therefore,  the  thermal 
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Figure  12.  Ra^fRa^^vs  Ma  in  the  case 
ofT2  =  4°C  (Seki  etal.  1977). 


Figure  13.  RaJ(R&c)Ma=0  vs  Ma  in  the  case 
o/Bi  =  0  (Seki  et  al.  1977). 


boundary  condition  of  the  free  surface  becomes  less  significant  as  T2  increases  and  the  effect  of  Ma 
on  Rac  is  reduced.  However,  the  effect  of  Ma  on  Rac  still  can  be  evaluated  even  at  T2  >  4°C.  This 
behavior  may  be  attributed  to  the  upward  penetration  of  the  free  convective  motion  in  the  unstable 
layer  exceeding  the  4°C  isotherm.  Figure  1 4  shows  the  comparison  of  experimental  and  analytical 
results,  in  which  Seki  et  al.  (1977)  evaluated  the  Revalue  based  on 

toe  =  gp  (rm-7'i)/iV  (17) 

for  T2  >  4°C  and 


Rac  =  sp  (r2-r,)W> a 
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(18) 


1977). 


for  T2  ^  4°C.  p,  v,  and  a  are  evaluated  at  2°C  in  eq  17  and  evaluated  at  the  arithmetic  mean  tem¬ 
perature  in  eq  18.  It  can  be  concluded  that  for  T2  <  8°C,  Rac  ~  500  agreed  extremely  well  with  the 
experimental  results  of  Sun  et  al.  (1969). 

In  a  melt  layer  between  vertical  plates 

The  only  analytical  and  experimental  work  dealing  with  the  onset  of  convection  of  a  vertical  water 
layer  formed  by  melting  ice  was  reported  by  Hassab  and  Sorour  (1982).  Their  problem  deals  with 
the  melting  of  an  ice  layer  that  is  confined  inside  a  vertical  slender  slot  as  shown  in  Figure  15. 
Initially  the  ice  is  at  its  melting  temperature  (i.e.,  T  =  0°C )  at  t  =  0,  then  the  left  side  is  suddenly 
increased  to  a  constant  temperature  Tj  >  0°C  and  maintained  at  that  value  throughout  the  entire 
experiment.  When  the  melt  layer  is  thin,  the  heat  is  transferred  by  conduction  (except  at  the  ends) 
and  consequently  a  laminar  parallel  flow  will  be  developed  as  a  result  of  ice  density  difference  in 
the  fluid.  As  the  melt  layer  thickness  gradually  increases  with  time,  the  initial  laminar  motion  breaks 
up,  and  secondary  flow  appears  in  the  form  of  either  stationary  horizontal  cells  or  traveling  waves. 

The  stability  problem  was  solved  by  the  application  of  the  Galerkin  method,  and  the  stability 
criteria  of  this  system  are  established  by  determining  its  eigenvalue,  Cn  =  Cr+  iCj  (where  C;  =  wave 
speed,  n  =  1, 2,....N).  That  is,  for  a  given  choice  of  the  system  parameter  T{,  there  are  at  least  two 
minimum  values  of  Gr  (Grashof  number)  with  respect  to  a  wave  number  (2 nfk,  where  X  is  the 
wavelength)  for  either  the  real  part  of  the  least  eigenvalue  Cr  =  0,  Ci  *  0,  in  the  case  of  traveling 
waves,  or  for  both  components  Cr = Cj = 0  in  the  case  of  stationary  cells.  Figure  1 6  shows  the  critical 
Grashof  number  along  with  the  modified  Grashof  number  Grc  as  a  function  of  7, .  The  values  of  Grc 
and  Grc  are  defined  respectively  as 
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Figure  15.  Schematic  of  the  problem  (Hassab 
and  Sorour  1982). 


20 


J _ L 

24  26 


'i  ( Hassab  and  Sorour  1982). 


and 


Grc  =  gy  i  T\  h3/v2 


(19) 


GrP  =  kshJL 


in  which  |30  is  evaluated  from  eq  4  as 


(20) 


Po  = 


_  J_ 

Po  dT 


T  =  0 


(2Y  l  —  2  Tm)  Tm . 


(21) 


Substituting  values  of  y(  =  0.793953x  1 0_5oC~2,  y2  =  -0.6559x  1 0~7oC“3  and  7m  =  4°C,  the  ratio  of  eq 
20  to  19  can  be  expressed  as 


Grc  =  8.3965  Grc/T\  .  (22) 

As  shown  in  Figure  16,  the  Grc  curves  approach  an  asymptotic  value  of  1 300  for  7,  >  26°C  and,  as  in 
horizontally  confined  melt  layers,  the  critical  Grashof  number  is  no  longer  a  constant  but  varies  with 
7 , .  However,  this  result  would  not  provide  insight  into  the  value  of  7,  at  which  the  instability  is 
enhanced  or  delayed.  This  uncertainty  is  due  to  the  dependence  of  Gr  upon  both  7,  and  the  melt  layer 
thickness  h(t).  It  is  therefore  more  appropriate  to  recast  the  stability  results  in  terms  of  critical  melt  layer 
thickness  hc  and  critical  melting  time  tc  which  are  derived  from  the  critical  Grashof  number.  The 
critical  melt  layer  thickness  and  critical  melting  time  are  defined  as 


hc 


v2  Grc 
Yl  gT  i 


(23) 


and 


h2 

tc  =  (24) 

4S2  a 

Figure  17  shows  the  variation  of  ftc  and  tc  as  a  function  of  7,.  The  graph  can  be  classified  into  three 
regions  as  follows: 

1.  In  the  region  when  7,  <  4°C,  the  density  profile  has  no  maximum  alue;  increasing  7, 
increases  the  density  difference  across  the  melt  layer,  accordingly  destabilizing  the  flow. 
When  7j  >  4°C,  the  melt  layer  has  a  maximum  density,  but  increasing  7,  has  a  smaller 
destabilizing  rate,  although  the  instability  sets  in  earlier  owing  to  the  higher  temperature 
difference  across  the  melt  layer. 

2.  In  the  region  for  7. 1  °C  <,  Tl  £  9.4°C,  the  variations  of  hc  and  tc  are  not  monotonous  for  in¬ 
creased  7 , .  This  trend  may  be  due  to  the  transition  in  the  wave  structure  from  two-column 
waves  for  7j  <  7. 1°C  to  three  column  waves  for  7. 1  °C  <7,  <  9.4°C  and  back  as  two-column 
waves  for  7,  >  9.4°C. 

3.  In  the  region  for  Tx  >  9.4°C,  increasing  7,  has  a  significant  destabilizing  effect  because  of 
the  strong  bending  of  the  density  profile  resulting  from  the  pronounced  increase  of  the 
density  difference. 

Figures  1 8a  and  b  show  the  dimensionless  mean  base  flow  velocity  across  the  melt  layer  thickness. 
Figure  1 8a  shows  the  U  distribution  in  regions  I  and  III.  In  region  I,  a  unicellular  motion  forms  so  that 


15 


01 

E 


o* 

c 


a> 

E 


o 

o 


o 


Figure  1 7.  Critical  melting  thickness,  hf  and  melting  time,  tc  as function  of step-wall  temperature  T  t 
(Hassab  and  Sorour  1982). 


the  fluid  near  the  hot  wall  descends  and  near  the  interface  ascends.  An  opposite  unicellular  motion 
occurs  in  region  III  because  the  water  near  the  interface  is  denser  than  that  near  the  hot  wall.  However, 
in  region  in  (Fig.  1 8b),  the  fluid  at  the  central  part  is  denser  than  that  next  to  the  boundaries.  Therefore, 
two  parallel  cells  form,  with  the  fluid  at  the  central  part  moving  downwards  and  that  adjacent  to  the 
Boundaries  moving  upwards.  Figure  19  shows  the  variation  of  average  volumetric  expansion  3  as  a 
function  of  7,  and  sheds  some  light  on  the  magnitude  of  U  in  regions  I  and  III.  When  7,  <  8.2°C,  P 
decreases  with  increases  of  7,;  for  7j  >  8.2°C,  P  increases  with  increasing  T j .  Since  the  shear  flow 
is  driven  by  the  buoyancy  effects  that  are  linearly  dependent  upon  p,  increasing  7,  will  suppress  the 
dimensionless  base  flow  U  in  region  I  while  expanding  it  in  region  III.  (See  Fig.  1 8a  to  compare  U 
curves  for  7,  =  4°  and  6°C,  and  for  7(  =  12°  and  20°C.) 

The  stability  problem  of  natural  convection  in  a  vertical  melt  layer  can  be  classified  in  terms  of  the 
boundary  temperature  7,  into  three  regions: 

1 .  For  7,  <  7. 1  °C,  stability  sets  in  as  two-column  waves  travel  opposite  in  the  vertical  direction. 
Although  heating  increases  Gr  it  is  actually  destabilizing  the  flow  with  a  decreasing  rate. 

2.  For  temperature  range  7. 1 0  <  7!  <  9.4°C,  instability  occurs  as  three-column  waves.  Minimum 
wavelength  and  minimum  wave  speed  occur  at  8°C  and  maximum  values  of  Grc  occur  at  7, 
=  7.3°  and  9.2°C,  respectively.  In  this  region  heating  also  destabilizes  the  flow. 

3.  For  7,  >  9.4°C,  the  instability  recurs  as  two  column  waves,  with  both  wave  speed  and 
wavelength  increasing  (see  Fig.  20  and  21)  as  7,  increases.  Again,  heating  has  a  destabilizing 
effect  on  the  flow. 
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a.  Regions  I  and  III. 


b.  Region  III. 


Figure  18.  Velocity  profiles  in  Regions  I  and  III  (Hassab  and  Sorour  1982). 
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EFFECT  OF  4°C  ON  HEAT  TRANSFER 
IN  PURE  WATER  AND  ICE-WATER  SYSTEMS 

In  a  vertical  enclosure 

The  vast  majority  of  the  natural  convection  studies  reported  so  far  pertain  to  fluids  whose  density 
decreases  or  increases  monotonically  with  the  temperature  at  a  constant  pressure.  However,  for  liquids 
such  as  water,  the  density-temperature  curve  exhibits  a  maximum  value.  Since  the  coefficient  of 
thermal  expansion  changes  sign  through  the  maximum,  the  use  of  linear  approximation  of  density  to 
model  water  is  inappropriate  in  the  range  of  temperature  that  corresponds  to  the  density  maximum. 
Lankford  and  Bejan  (1986)  were  the  first  to  investigate  the  natural  convection  in  a  vertical  enclosure 
filled  with  water  near  4°C.  Two  series  of  experimental  runs  were  reported:  one  set  for  rH  (hot-wall) 
and  Tc  (cold-wall)  temperatures  well  above  Tm (temperature  corresponding  to  maximum  density),  and 
one  set  for  7"H  and  Tc  embracing  Tm. 

Figures  22  and  23  show  typical  wall  temperature  distributions  as  a  function  of  dimensionless  height 
y/H  and  their  corresponding  flow  patterns.  When  the  temperature  everywhere  in  the  enclosure  is  well 
above  4°C  the  wall  temperature  increases  almost  linearly  with  height  and,  at  each  level,  the  wall-to- 
wall  temperature  difference  is  practically  constant  (Fig.  22).  In  experiments  where  the  entire  enclosure 
is  below  4°C,  the  water  density  consistently  increases  upon  heating.  The  wall  temperature  distribution 
in  such  a  case  is  similar  to  that  shown  in  Figure  22  except  that  the  temperature  of  both  walls  decreases 
with  height. 

When  the  enclosure  temperature  range  includes  4°C  and  when  the  average  enclosure  temperature 
is  slightly  less  than  4°C,  the  sinking  core  jet  emerges  from  near  the  top  of  the  warm  wall,  and  along 
the  cold  wall  the  water  jet  rises  all  the  way  to  the  top  as  it  is  cooled  to  a  temperature  close  to  0°C.  The 
same  jet  turns  around  at  the  top  of  the  enclosure  and  descends  along  the  warm  wall,  because  its  density 
increases  as  its  temperature  increases  from  0°  to  4°C. 
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Figure  22.  Example  of  wall  tem¬ 
perature  distribution  withTHand 
T  consistently  greater  than  Tw 
(Lankford  and  Bejan  1986). 


Figure  23.  Example  of  wall 
temperature  distribution  when 
Th  and  Tc  embrace  Tm  (Lank¬ 
ford  and  Bejan  1986). 


When  the  average  enclosure  temperature  is  slightly  above  4°C  the  sinking  core  jet  emerges  from  a 
point  near  the  top  of  the  cold  wall.  Along  the  warm  wall,  which  is  swept  by  a  jet  whose  density 
decreases  monotonically  with  temperature,  the  temperature  distribution  has  the  same  features  as  those 
in  an  enclosure  heated  with  uniform  heat  flux  from  the  side  and  filled  with  fluid  whose  density 
conforms  to  the  linear  density-temperature.  Along  the  cold  wall,  however,  the  collision  of  the  two 
counterflowing  jets  is  felt  as  a  sharp  change  in  vertical  temperature  gradient  near  the  level  of  collision 
(see  Fig.  23). 

Figure  24  shows  the  experimental  results  for  the  case  of  an  average  water  temperature  well  above 
4°C  along  with  empirical  correlations  recommended  for  the  cases  of  isothermal  walls  and  constant- 
heat-flux-walls.  In  this  case,  the  heat  transfer  data  can  be  well  represented  by 

Nu  =  0.203  Ran29  (25) 

in  which 

Nu  =  - -g-—  (26) 

Wk  [Th  -  Tc) 

and 

(27) 

va 

where  Q  =  the  overall  heat  transfer  rate, 

W  =  the  enclosure  width, 

Th  =  the  height, 

T c  =  height-averaged  wall  temperatures. 


Figure  24.  Correlation  of  heat  transfer  results  in  water  at  temperatures  above  4  °C  (Lankford 
andBejan  1986). 
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Figure  25.  The  breakdown  of  the  Nua/idRa^  correlation  ( eq25)asTm  is  approached  ( Lankf ord 
and  Bejan  1986). 


The  physical  properties  k,  a,  p,  and  v  were  evaluated  at  the  average  enclosure  temperature  defined  as 

Tavg  =  £  {Th+Tc)  (28) 

The  experimental  results  are  found  to  agree  well  with  the  isothermal  wall  correlation  reported  by  Bejan 
(1979),  i.e., 

Nu  =  0.364  Rail 4  .  (29) 

However,  the  slope  of  the  empirical  correlation  of  eq  25  is  nearly  identical  to  that  of  eq  30  reported 
by  Kimura  and  Bejan  (1984)  for  the  case  of  constant-heat-flux  wall,  i.e., 

Nu  -  0.25  (If.)'"  Raff .  <»» 

Figure  25,  a  plot  of  Nu/RaH0  2S  vs  Tc ,  shows  clearly  how  these  correlation  methods  deteriorate 
as  temperatures  approach  4°C.  The  maximum  density  flow  pattern  inhibits  the  transfer  of  heat, 
causing  a  decrease  in  Nusselt  number  that  was  not  accounted  for  in  the  development  of  eq  25.  Even 
in  the  range  of  temperatures  immediately  above  4°C,  where  the  density  maximum  has  not  neces- 
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Figure  26.  Schematic  of  the  boundary  layer  regime 
in  a  vertical  space  filled  with  cold  water  at  a  mean 
temperature  of4°C  (Lankford  and  Bejan  1986). 
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sarily  been  crossed,  the  data  still  tend  to  deviate  from  the  constant  value  assumed  by  eq  25  (see  the 
dotted  line  in  Fig.  25).  This  is  because  the  coefficient  of  thermal  expansion  p  used  in  the  linear 
density-temperature  approximation  (i  .e.,  p = p0  [  1  -P(T-ro)]  where  po  and  Tq  are  the  reference  density 
and  temperature,  respectively),  changes  so  rapidly  near4°C  that  averaging  it  in  the  normal  fashion 
leads  to  significant  errors. 

To  overcome  this  difficulty,  Lankford  and  Bejan  ( 1 986)  adopted  the  technique  of  scale  analysis  by 
considering  the  boundary  layer  regime  for  natural  convection  in  near  4°C  water.  Figure  26  shows  the 
case  where  the  side-wall  temperatures  occur  symmetrically  above  and  below  4°C.  It  shows  that  the 
maximum  density  sinking  jet  originates  from  near  the  top  of  the  enclosure.  The  convection  pattern 
consists  of  two  boundary  jets  rising  along  the  differentially  heated  side  walls,  and  a  maximum  density 
sinking  core  with  Tm  =  4°C.  Thermal  boundary  thickness  8H  and  Sc  can  be  scaled  from  the  classical 
case  for  fluid  with  Pr>  1  (Bejan  1984),  i.e., 

5h  -  P/P)  1/4  (31) 

V  va  / 

Expressing  Ap/p  =  y(T  -Tm)2  (from  eq  2),  eq.  31  becomes 

6„  .  (32) 

Similarly  for  the  cold-side  boundary  layer  we  have 

va 

The  quantity  in  the  square  brackets  in  eq  32  and  33  is  the  Rayleigh  number  for  fluids  with  parabolic 
density  variation  near  the  density  maximum. 
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Figure  27.  Correlation  of  heat  transfer  results  in  water  at  temperature 
near  4  °C  ( Lankford  and  Bejan  1 986 ). 

The  scale  of  overall  heat  transfer  rate  from  th  to  Tc  can  be  written  as 

q  ~  Ih.-  Jc. 

Rh+Rc 

where  /?H  ~  &H/kHW  and  Rc  =  &JkHW.  From  eq  32  and  33,  eq  34  can  be  expressed  as 

Nu  =  - - 

RaiH  +  C2Raic 

where 


(34) 


(35) 


Ray h 


gyH3  {TH-Tm)2 

~V  * 


Rfl)C  - 


g7tf3(7'm-7-c)2 

va 


(36) 


Figure  27  shows  an  excellent  correlation  of  the  heat  transfer  results  of  water  at  near  4°C  with  eq 
35  derived  from  scale  analysis.  The  solid  line  drawn  through  the  experimental  points  is  a  plot  of  eq 
35  with  assumed  values  of  C,  =  0.3 1  and  C2 = 0.5.  These  constants  were  chosen  such  that  the  standard 
deviation  between  eq  35  and  the  measurements  was  minimized  when  summed  over  all  the 
experimental  points. 

In  a  vertical  annulus 

Lin  and  Nansteel  (1987b)  were  the  first  to  investigate  systematically  the  effect  of  curvature  and 
aspect  ratio  on  the  flow  structure  and  heat  transfer  in  a  vertical  annulus.  Figure  28  shows  the  vertical 
annulus,  coordinate  system  and  thermal  boundary  conditions.  Lin  and  Nansteel  ( 1 987b)  assumed  that 
the  flow  in  the  gap  is  laminar  and  axisymmetric,  and  all  fluid  properties  are  constant  except  for  density 
in  the  buoyancy  term  of  the  vertical  momentum  balance.  They  introduced  dimensionless  variables 
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Figure  28.  Vertical  annulus,  coordinate  system  and  thermal 
boundary  conditions  (Lin  and  Nansteel  1987b). 


R  =^— 5-,  Z  =X,  (/  =  «£,  V  =?£- 

D  D  vc  vc 


x  =^sl,  p  a  P°2  ,  e  ~—~J P- 

D2  PcV2  T’h-T'c 

and  represented  the  density-temperature  relation  as 

P  -  Ph.[i-a,|r-rJ«].  (37) 

The  governing  equations  for  conservation  of  mass,  momentum  and  energy  in  the  annulus  are  given  as 
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DQ  _  l  d2e  +  (K-  0  96  +  9^ 

Dt  Pr[dR  2  (jf-l)*  +  l  9 R  dZ2\' 


where  U  =  V  =  0  on  the  boundary 
6  (0,  Z,)  =  1,  0(1,  Z)  =  0 

(fl,0)  =  (R,A)  =  0. 

9Z  9Z 

K  =  the  ratio  of  rjr.^,  which  characterizes  the  degree  of  curvature, 

Pr  =  the  Prandtl  number,  v  Jac 
A  =  the  aspect  ratio  H/D. 

The  Rayleigh  number  is  defined  as 

Ra'  =  a*  (^h  ~  ^2) 

Pc  vc  ac 

and  the  density  distribution  parameter  as 

r'  -  Tm~  Tc  (43) 

Th-Tc 

which  fixes  the  orientation  of  the  maximum  density  temperature  with  respect  to  the  wall  temperatures 
Th  and  r  .  For  example,  R'=  1/2  corresponds  to  the  circumstance  in  which  the  hot  and  cold  wall 
temperatures  perfectly  straddle  the  maximum  density  temperature  for0<i?'<  1,  (i.e.,  Th>  Tm),  and 
the  water  density  in  the  annulus  increases  with  increasing  temperature  to  a  maximum  P  =  Pm  at  T  = 
Tm  and  then  decreases  with  any  further  increase  in  temperature.  On  the  other  hand  R'=  0  corresponds 
to  the  case  in  which  the  cold  wall  is  at  temperature  Tm  so  that  fluid  density  decreases  monotonically 
with  temperature  everywhere  in  the  annulus.  For  R  '=  0,  flow  in  the  annulus  is  similar  to  the  flow  of 
a  Boussinesq  fluid,  but  only  in  a  qualitative  sense  since  density  is  a  nonlinear  function  of  temperature. 
When  R'=l,Th  =Tmthe  density  increases  monotonically  with  temperature  everywhere  in  the  annulus. 
Eliminating  P  in  eq  30  and  40  and  introducing  the  dimensionless  vorticity  t, 

e  _  9V_  _  dU_ 

5  9 R  9 Z 

and  stream  function  V 


u=  (g-0  a*. 

(if  -  1 )  /?  +  1  9Z  ’ 


{K- 1)  3\j/ 

(tf-l)/?  +  l  dR 


we  have 
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with 


t  = _ 1  (d2  V  +  _  (*~0  dyA 

"  {K-l)R  +  l  \a/?2  dz2  (K-\)R  +  l  dR  I 

The  Nusselt  numbers  at  the  inner  and  outer  walls  can  be  expressed  as 


Nui  = 


k(Th-Tc) 


* A 

o 


^.(o,  Z)dZ 
d  R 


and 


Nu0  = 


k(Th-Tc) 


(1  ,Z)dZ. 


(46) 


(47) 


Nu.  is  related  to  Nu  by 

1  O  J 

Nu;  =  KNu0.  (48) 

Numerical  results  were  obtained  by  solving  eq  41  and 45  with  the  finite  difference  method.  Figure 
29  shows  the  effect  of  R'  on  flow  structure  and  temperature  distribution  for  the  case  of  A=1  (aspect 
ratio)  andAT=l  (no  curvature  effect)  and  Ra'-  105.  In  Figure  29a,  R  '=  0.4,  the  two  counterrotating 
cells  are  separated  roughly  by  the  maximum  density  isotherm  0  =  R' =  0.4.  In  Figure  29b,  the 
relatively  light  fluid  rises  adjacent  to  the  heated  and  cooled  walls,  while  dense  fluid  p  =  pm,  T  ~  Tm 
falls  near  the  enclosure  vertical  midplane  in  a  symmetric  dual-cell  pattern.  As  /?' increases,  the 
counterclockwise  rotating  cell  has  grown  substantially  in  size  and  strength  at  the  expense  of  the 
clockwise  rotating  cell.  In  the  case  of  R'=  1,  the  density  increases  with  temperature  everywhere, 
resulting  in  a  single  counterclockwise  vortex.  This  flow  pattern  is  qualitatively  similar  to  the 
convection  of  a  Boussinesq  fluid  except  for  the  direction  of  circulation. 


a.  R'  =  0.4. 

Figure  29.  Stream  Junction  x  10 f)  and  temperature  (0  =  0(0.1  )1)  contours  for  A  =  /,  K  = 

1,  and  Ra'=  105  (Lin  and  Nansteel  1987b). 
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Figure  29  (cont’d).  Stream  function  (VfxlO3)  and  temperature  (6  =  0(0. 1  )1)  contours  for 
A  =  /,  K  =  /,  and  Ra'=  JO5  (Lin  and Nansteel  1987b). 
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Figure  30  shows  the  effect  of  K  on  flow  structure  and  temperature  contours  for  A  =l,Ra'=  105,  and 
R  '=  0.4.  As  indicated  in  eq  46,  the  average  heat  flux  at  the  inner  walls  is  AT  times  larger  than  the  average 
heat  flux  at  the  outer  wall  under  steady  conditions.  As  a  result,  with  increasing  K,  isotherms  will  become 
increasingly  crowded  near  the  inner  cylinder.  As  the  maximum  density  water  shifts  toward  the  inner 
cylinder,  so  does  the  line  of  demarcation  between  the  outer  counterclockwise-rotating  cell  and  the 
inner  clockwise-rotating  cell.  Thus  for  fixed  R ",  the  outer  cell  strengthens  and  the  inner  cell  weakens 
with  increasing  curvature  of  the  annulus  (i.e.,  as  K  increases). 

Figure31  shows  the  effect  of  A  on  y  and  0  for  the  case  ofK-2,R'=  1/2  and  Ra'-  104.  Clearly, 
increasing  A  tends  to  increase  the  convective  intensity  of  the  outer  cell.  This  is  because  the  quenching 
effect  of  the  unheated  horizontal  surfaces  lessens  as  A  is  increased  and  the  temperature  field  becomes 
dominated  by  conduction  over  much  of  the  vertical  span  of  the  annulus  when  A  >  4.  Only  near  the  top 
of  the  annulus,  where  the  counterclockwise  outer  cell  sweeps  cool  fluid  from  the  outer  wall  onto  the 
heated  inner  wall,  does  the  contribution  of  convection  become  significant.  This  phenomenon  is  not 


Figure  31 .  Effect  of  A  on  distribution  of  stream  function  (yxlO3)  temperature  contour  ( 8 
=  0(0.1  )1)  for  R'  =  0.5,  K  =  2  and  Ra'  =  10*  (Lin  andNansteel  1987b). 
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Figure  31  (cont'd). 


observed  near  the  bottom  of  the  annulus  space  because  of  the  existence  of  a  dual-cell  structure,  which 
prevents  the  contact  between  the  inner  cell  and  outer,  cooled  wall. 

Heat  transfer 

Figure  32  shows  the  variation  of  Nu{  with  K(K  =  rjr)  for  the  case  R'=  0.4,  A  =  1  and  a  range  of 
103  <*Ra'<,  105.  In  this  particular  condition,  there  is  a  transition  from  inner  to  outer  cell  dominance 
as  curvature  (i  .e.,  as  K )  increases.  This  occurs  in  the  range  of  2  <  K  <  5  for  1 03  £  Ra'£  1 05.  The  smooth 
transition  takes  place  from  a  case  where  the  inner  cell  wets  both  walls  to  one  in  which  the  outer  cell 
wets  both  walls.  Over  much  of  this  range  of  K,  the  cells  are  of  roughly  equal  strength  and  the  inner  cell 
is  effectively  insulated  from  the  outer  wall  by  the  outer  cell  and  vice  versa.  This  in  turn  leads  to  a 
diminished  rate  of  convective  heat  transfer  from  the  inner  to  the  outer  wall.  The  result  is  a  minimum 
heat  transfer  rate  near  K  -  3.5,  but  only  a  minimum  rate  of  increase  of  Nu-  with  K,  since  conduction 
heat  transfer  increases  quite  rapidly  with  increasing  K.  The  curve  for  Ra'=  103  shows  no  point  of 
inflection  because  the  contribution  of  convection  to  Nui  is  rather  small.  The  variation  of  (/V«j)condcan 
be  related  to  K  by 


(A^cond  =  VT-  (49) 

In  a 

The  value  of  (/V«.)cond  is,  in  general,  not  in  unity  as  it  is  in  the  case  of  K  =  1  (rectangular  enclosure). 
The  effect  of  K  and  Ra  'on  Nu.f  is  easily  observed  by  plotting  Nu.J{Nu.^coni  vs  AT  as  it  is  shown  in  Figure 
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K 


Figure  32.  Variation  o/Nu-  with 
KforR'  =  0.4,  A  =  1  using  Ra'as 
a  parameter  (Lin  and  Nansteel 
1987b). 


g  -j  8  9  10  Figure  33.  Variation  of  Ratio  of 

Nui/(Nui)condforR'  =  0.4,A  =  I 
K  (Lin  and  Nansteel  1987b). 


33.  This  figure  shows  the  convective  contribution  to  the  decrease  of  heat  transfer  with  increasing  K 
for  K  <,  3.5  and  then  the  increase  for  larger  K  (behavior  that  is  a  direct  consequence  of  the  transition 
from  inner  to  outer  cell  dominance).  For  large  values  of  K,  the  ratio  again  decreases  gradually  since 
conduction  heat  transfer  grows  without  bound  while  convective  heat  transfer  is  bounded.  Conse¬ 
quently,  all  the  curves  will  approach  a  value  of  unity  asymptotically  as  K  — ►  «>. 

Between  horizontal  concentric  cylinders 

Sekietal.  (1975)  reported  their  pioneer  work  on  natural  convection  heat  transfer  between  horizontal 
concentric  cylinders  of  water  near  its  maximum  density.  In  their  study  the  ratio  of  outer  to  inner 
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diameter  (djd.^  ranged  from  1.18  to  6.39  with  values  of  19.0, 38.0, 55.4, 69.6  and  99.7  mm  mdd0 
of  65.5  and  121 .5  mm.,  and  the  gap  width  L  [=  (d0-d 4)/2]  was  varied  from  5.0  to  5 1 .2  mm.  The  surface 
temperature  of  the  inner  cylinder  was  kept  atO°C,  whereas  the  outer  cylinder  surface  was  varied  from 
1°  to  15°C.  For  each  geometrical  configuration,  two  sets  of  data  were  obtained  (i.e.,  one  set  for  T0  < 
4°,  and  the  other  set  for  T0  >  4°C). 

Flow  patterns 

Figures  34  a,  b,  c,  and  d  show  a  series  of  photographs  and  schematic  views  of  flow  patterns  with  the 
sa med0  =  121.5  mm  but  with  increasing  d{  values  of  38.0  to  55.4,  69.6  and  99.7  mm  (i.e.,  with  djdt 


b.6i  =  55.4mm,T0  =  5°C. 


Figure  34.  Photographs  and  schematic  views  of  flow  patterns  for  dc  =  121.5  mm  (Seki  et  al.  1975). 
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Figure  34  (cont’t).  Photographs  and  schematic  views  of flow  patterns for  d0  =  121.5  mm  (Seki  etal. 
1975). 


ranging  from  3.2  to  1 .22  and  L  varying  from  41 .7  to  10.9).  The  temperature  of  T0  increased  from  4°C 
in  Figure  34a,  to  T0 = 5°C  in  Figure  34b,  and  T0  =  6°  in  Figure  34c;  however,  T0 = 4°C  in  Figure  34d, 
as  in  Figure  34a.  In  Figure  34a,  the  flow  patterns  are  stable,  with  the  eddy  moving  upward  along  the 
inner  cylinder  and  downward  along  the  outer  cylinder.  The  water  near  the  upper  part  of  the  vertical 
axis  flows  at  higher  speed  than  that  on  the  other  part  in  the  gap.  The  center  of  the  eddy  is  clearly 
observed  in  the  upper  part  of  the  annuli,  while  in  the  lower  part  a  secondary  counter-eddy  can  be  seen. 
For  T0 = 5°C  and  djdx  =  219  (Fig.  34b),  the  secondary  counter-eddy  in  the  lower  part  can  be  seen  more 
clearly .  When  T0  ranges  from  6°  to  7°C,  the  influence  of  density  inversion  at  4*0  becomes  more  evident. 
The  eddy  becomes  large  and  finally  two  counter-eddies  of  almost  equal  strength  exist  in  the  gap  (Fig. 
34c).  Flow  patterns  for  L  =  10.9  mm  and  djdx  =  1.22  with  T0  =  4°C  (Fig.  34d)  are  considerably 
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Figure  35.  Photograph  and  sche¬ 
matic  view for  dQ = 65.5  mm,  d; = 19.0 
mm  andT0  =  8°C  (Seki  et  at.  1975). 


different  from  those  in  Figure  34a-c.  A  very  small  eddy  that  flows  downward  along  the  vertical  axis 
is  observed  in  the  vicinity  of  the  top  annuli  for  all  T0.  This  is  believed  to  be  due  to  the  effect  of  L 
and  the  cylinder  curvature.  Figure  35  shows  the  case  of  dQ  =  65.5  mm,  dt  =  19.00  mm  (dQ/di  =  3.44), 
L  =  23.2  and  T0  =  8°C.  Here  the  eddy  is  flowing  upward  along  the  outer  cylinder  and  gradually 
extending  its  size,  while  the  one  moving  upward  along  the  inner  cylinder  becomes  smaller  and 
finally  disappears  as  T0  increases. 

Temperature  distribution  and  local  Nusselt  number 
Figures  36a-38a  show  the  effect  of  T0on  the  dimensionless  temperature  distribution  0  as  a  function 
of  the  dimensionless  radius  R.  The  corresponding  local  Nusselt  number  N«loc  as  a  function  of  angular 
position  (measured  from  the  upper  vertical  axis)  is  shown  in  Figures  36b-38b.  The  local  /Vn(oc  is  defined 
as 


Nu  ioc 


KAT 


(50) 


where  q  is  a  local  heat  flux.  Different  means  were  used  to  evaluate  this  quantity  to  determine  the  local 
Nusselt  number  for  the  inner  cylinder  (Auloc  j)  and  outer  cylinder  (Au)oc  0).  In  Figure  36a,  a  tem¬ 
perature  inversion  is  observed  in  the  range  from  0  to  90°,  which  is  reduced  as  angular  position  in¬ 
creases.  In  the  160-180°  range,  little  convective  flow  occurs.  Figure  36b  shows  a  minimum  Am  at  0° 
on  the  inner  cylinder  and  180°  on  the  outer  cylinder.  This  phenomenon  fits  well  with  the  observed 
temperature  distribution  in  Figure  36a.  Figure  37b  shows  the  case  of  T0  =  6°C,  where  two  standing 
eddies  of  almost  the  same  size  coexist.  The  An|oc  0  has  its  lowest  value  at  around  90°  angular  position, 
while  the  increase  of  Am)oc  ( with  increasing  angular  position  is  minimal  at  OP  and  highest  at  1 80°.  Figure 
38a  is  the  case  when  T0  =  8°C;  here  the  temperature  at  0°  angular  position  is  the  highest  one  in  the  gap 
between  the  cylinders  and  it  decreases  as  the  angular  position  increases.  In  the  range  of  angular  position 
from  0°  to  30°,  temperature  is  much  higher  than  those  in  either  angular  positions  and  it  implies  that 
there  is  a  stagnant  region  where  the  conduction  heat  transfer  is  dominant.  The  local  Nusselt  number 
for  the  outer  cylinder  decreases  as  the  angular  position  decreases  and  reaches  a  minimum  value  at  0°. 
On  the  other  hand,  the  Nuloc  i  values  decreases  gradually  from  90°  to  1 80°,  while  between  90°  to  0° 
it  increases  significantly  and  reaches  a  maximum  at  0°  (See  Fig.  38b).  Figure  39  shows  the  variations 
of  average  Nusselt  number  Nu  (based  on  inner  diameter)  with  AT (=  T0).  The  Nu  is  defined  as 
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a.  Temperature  profiles. 


b.  Local  Nusselt  numbers. 

Figure  36.  Temperature  profiles  and  local  Nusselt  numbers  for  AJ 
d,  =  3.20  and T 0  =  4°C  (Seki  et  al.  1975). 


Nu  =  —Oik. —  (51) 

ndi  ATk 

where  Q{  is  total  heat  transfer  rate  per  unit  length  of  the  cylinder.  Figure  39  shows  that  Nu  does  not 
increase  monotonically  with  increasing  A T  (as  is  common  in  fluids)  when  a  single  large  eddy  occupies 
the  major  portion  of  the  gap,  as  in  the  case  of  T0  £  4°C  or  T0  £  9°C.  The  peak  value  of  Nu_  is  found  at 
~4°C  and  the  minimum  values  appear  between  approximately  6°  and  7°C.  Values  of  Nu,  including 
minimum  Nu  values,  generally  increase  with  increasing  gap  width  L. 
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180* 

b.  Local  Nusselt  numbers. 

Figure  37.  Temperature  profiles  and  local  Nusselt  numbers  for  d = 
1.75,  and! 0=  6°C  (Seki  et  al.  1975). 


Nguyen  et  al.  ( 1 982)  were  the  first  to  conduct  an  analytical  study  on  the  natural  convection  of  cold 
water  between  two  horizontal  concentric  cylinders  with  constant  surface  temperatures  at  low 
Rayleigh  numbers.  The  governing  equations  were  solved  by  the  perturbation  method  and  the 
solutions  of  dimensionless  temperature  6'  and  stream  function  were  expressed  as  a  power  series 
of  the  nonlinear  Rayleigh  number  Ran  as 
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a.  Temperature  profiles. 


b.  Local  Nusselt  numbers. 


Figure  38.  Temperature  profiles  and  local  Nusselt  numbers  for 
cyd-  =  3.44  and  T a  =  8°C  (Seki  et  al.  1975). 


AT  (s T.) 


Figure 39.  Variation  of  Nu  with  AT  (or  TQ)  (Seki  etal.  1975). 


where  0'  is  defined  as  (T~To)f(T  -TQ),  and  Ran  is  defined  as 
Ran  =  g^(T.-To)^a 

in  which  P2  is  the  volumetric  coefficient  of  expansion  of  the  following  density-temperature  relation 

p  =  Pr[i  -  Pi  (r- r,)  -  fo(T-Ttf]  (54) 

where  7^ and  pfare  the  reference  temperature  and  density  respectively.  p(  is  related  to  P2  by 

PrWi 

with  p2  =  8x1  CT60^2  and  Tm  =  3.98°C. 

Although  the  disadvantage  of  eq  54  is  its  limited  use  for  the  range  of  temperatures  from  0°  to  8°C, 
an  advantage  is  that  this  equation  can  predict  all  the  essential  features  associated  with  an  inversion  of 
density  by  introducing  just  one  more  parameter  than  the  equation  for  a  classical  normal  fluid.  Figures 
40a-40c  show  the  effect  of  inversion  parameter  y  defined  as  P(/P 2(^r^o) or  (Tm-T0)/(T-T0)  on 
flow  patterns  and  isotherms  in  which  the  increments  between  adjacent  isotherms  and  streamlines  are, 
respectively, 

A0  =0.2 

5 

and 
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a.Y  =  0  and  Vmax  =  4.6. 


b.  y'=  -1.  w  =  -0.045 and  y  .  =-0.71. 

1  Tmax  Tmin 

Figure  40.  Streamlines  and  isotherms  for  Ra„ 


=  2000,  R"  =  2  (after  Nguyen  etal.  1982). 


AV=dVmaxl  +  l¥minl)/5, 

¥max  ¥„,„  being  the  values  of  the  stream  functions  at  the  centers  of  the  clockwise  and  counter¬ 
clockwise  vortices,  respectively. 

Figure  40a  shows  the  case  when  the  density  maximum  is  situated  at  the  outer  cylinder  (i.e.,  y ' = o). 
The  flow  consisted  of  two  symmetrical  counterrotating  vortices  with  a  downward  motion  near  the 
outer  cylinder.  The  resulting  flow  pattern  is  of  a  tadpole  shape  lying  in  the  upper  part  of  the  annulus 
where  the  fluid  motion  is  strongest.  The  maximum  heat  transfer  occurs,  however,  both  at  the  top  of 
the  outer  cylinder  and  at  the  bottom  of  the  inner  cylinder.  Since  there  is  no  density  inversion  in  this 
case,  the  flow  is  similar  to  that  observed  in  an  ordinaiy  fluid. 


Figure  40b  is  the  case  of  complete  inversion  with  y '  =  -1 .  The  heavy  dashed  line  represents  the 
4°C  isotherm  corresponding  to  the  region  of  maximum  density.  In  the  neighborhood  of  this  region, 
the  fluid  moves  downward,  while  near  both  the  inner  and  outer  cylinders  the  fluid  moves  upward. 
The  effect  of  the  inversion  is  the  sharp  decrease  in  heat  transfer  as  indicated  by  the  almost  concentric 
isotherms  in  contrast  with  the  rather  distorted  patterns  shown  in  Figure  40a.  Figure  40c  is  the  case 
of  Y  =  -2,  (i.e.,  when  the  maximum  density  is  situated  at  the  inner  cylinder)  where  the  flow  and 
isotherm  patterns  are  just  opposite  to  those  observed  in  Figure  40a. 

Figure41ashows  the  angularvelocity  profiles  at0  =  90°for/?an=  8000, f?"=2withy'asa  function 
of  (r/r j  - 1 ).  Figure  4 1  b  provides  the  variation  of  angular  velocity  as  a  function  of  angular  position. 

The  effect  of  /?anon  the  flow  and  isotherms  pattern  can  be  seen  by  comparing  Figures  42a  and  40b. 
By  increasing  Ran,  the  inner  cell  is  pushed  down  while  the  outerone  moves  up  and  gradually  changes 
from  the  well-  known  kidney  shape  to  that  of  a  tadpole.  The  heat  convection  is  enhanced  in  such  a 
way  that  the  isotherms  are  squeezed  to  the  top  of  the  inner  cylinder  and  to  the  bottom  of  the  outer 
one.  The  effects  of  changing  radius  ratio  [/?' = (rQ/r)]  on  the  inversion  phenomenon  can  be  seen  from 
Figures  42a  and  42c  or  42b  and  42d. 

An  average  or  overall  Nusselt  number  can  be  defined  as 

Nu  =  1  +  e  Roq  (55) 
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R"= /.2Y=-7.Vmax=0.00S7,andVmin=-0.0/3.  d.  R"  =  2.6,  y'  =  0.557,  vmax  =  3.2  andx ^ 


Figure  42  (coni'  d). 


where  RaQ  is  the  Rayleigh  number  based  on  gap  width  (rQ-r),  i.e., 

5«G  =  gP2(ro-ri)3(ri-r)2/va 

where  P2  =  3x  1 0-6  °C"2  and  e  is  a  function  of  R"  and  y'.  Equation  55  explicitly  shows  that  the  overall 
Nusselt  number  is  made  of  two  terms,  the  first  term  representing  heat  transfer  due  to  pure  conduction 
and  the  second  arising  from  heat  transfer  by  convection.  Figure  43  shows  the  effects  of  density 
inversion,  for  R"  values  ranging  from  R"  =  1.25  to  10,  on  the  overall  convective  heat  transfer 
coefficient  e  =  ( Nu-\)IRaG 2.  All  three  curves  present  a  very  sharp  minimum  in  the  vicinity  of  y '  = 
-1,  indicating  a  deep  cut  of  the  heat  transfer  rate  due  to  the  complete  inversion  of  fluid  density.  For 
example,  the  minimum  value  of  e  at  y'  =  -1  is  about  103  times  smaller  than  the  value  of  e  at  y '  = 
-2.  This  steep  drop  is  essentially  due  to  the  presence  of  two  counterrotating  vortices  that  arise  from 
the  inversion  of  the  fluid  density. 

Vasseur  et  al.  ( 1 983)  extended  Nguyen  et  al.’s  ( 1 982)  and  Seki  et  al.’s  (1975)  experimental  work 
by  conducting  an  analytical  study  of  the  same  problem.  Vasseur  et  al.  (1983)  assumed  symmetry 
along  the  vertical  plane  through  the  axis,  neglecting  the  viscous  dissipation  and  compressibility 
effects.  All  fluid  properties  except  the  water  density  were  taken  to  be  constant  and  evaluated  at  the 
arithmetic  mean  temperature  of  the  two  cylinders.  The  governing  equations  for  the  present  problem, 
using  the  Oberback-Boussinesq  approximation,  are 


dn  +  J_  + 

dx  R"  dRm 


=  PrV2a.  +  A 


sin<|> 


9Ag_  + 
dR"' 


cos<}>  _  dAp 
R'"  3<t>  . 


(56) 
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Figure  43.  Coefficient  of  convective  heat  transfer  £  as  function  of  y '  and  R"  (Nguyen  et  al. 
1982). 


U  =  -L  V  =  - 

R’"  30  ’  3 R"'  ' 

The  initial  and  boundary  conditions  are 

U=V  =  =  0;  0  =  0,  at  t  =  0 


U=V  =  \/  =  0;  0  =  0,  0  <  0  <  K,  R'" 


1 

R"  -  I 


U=V  =  y  =  0;  0=1,  0<  0  <  71,  R"' 


—  =  V  =  y  =  Q  =  0;  —  =  0  at  0=0,  n  (symmetry  lines) 

30  30 

where  fi  =  dimensionless  vorticity,  (aL2/a 
©  =  the  vorticity 
L  =  the  gap  width  (r-r.) 
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dimensionless  radial  coordinate,  r/L 

dimensionless  radial  and  angular  velocity,  uL/a  and  vL/a 

radial  and  angular  velocity 

the  thermal  diffusivity 

the  angular  coordinate, 

defined  as  [p  -  p(0)]/p 

the  water  density 

an  average  water  density. 

dimensionless  temperature,  ( T-T-^/AT  (A T  =  7^-jT) 
dimensionless  stream  function  (=  y/a) 

the  stream  function 
radius  ratio  (rQ/r). 


This  v.  ^rk  covers  a  temperature  range  from  0°  to  20°C  with  the  density-temperature  relation  rep¬ 
resented  by 


^  =  1  +  pjr  +  T2+  p/3  +  p474  (61) 

where  po  =  0.9998396  (g/cm3)  is  water  density  at  0°C,  with  P'  values  given  by 
p'°  =  -0.6789452x10-4(1/°C) 

P'2  =  0.907294338x  1 0~5(  1/°C2) 
p'3  =  0.964568 1 25x  1 0~7(  1  /°C3) 
p'4  =  0.873702983xl0~9(l/°C4). 

A  nonlinear  Rayleigh  number  Ran  and  an  inversion  parameter  y '  were  used  to  characterize  their 
computed  results: 


Ra„  = 


—A _  h  A T2 

(R"-lf  Pr 


(62) 


and 


y'  =  2  lTm~  la)  (63) 

'  A  T  I 

where  A  =  a  size  parameter  defined  as  gL3/ a2 
P2  =  a  coefficient  defined  in  eq  6 1 
AT  =T-T. 

O  1 

and  y '  is  related  to  R'  defined  in  eq  43  by  y '  =  -2 R'. 

In  their  study,  the  value  of  y'  varies  between  -2  S  y '  ^  0.  When  y'  =  -1 ,  Tm  is  half  way  between 
TQ  and  Ti  and  the  densities  on  the  inner  and  outer  cylinder  are  the  same.  When  y'  <-2  and  y'  >  0, 
no  inversion  occurs.  When  ly'l »  1  the  fluid  behavior  tends  to  that  of  a  common  fluid  with  a  linear 
equation  of  state.  Figures  44a-44d  show  the  variation  of  the  isotherm  and  streamlines  as  a  function 
of  Ran,  R"  and  y'.  In  Figure  44a,  y ' = 0  (i.e.,  ro-4°C),  and  the  fluid  near  the  outer  cylinder  is  heavier 
and  therefore  is  moving  downward  while  the  relatively  higher  fluid  near  the  inner  cylinder  is  moving 
upward  (this  is  nearly  equivalent  to  Fig.  34a  and  36a).  The  distortion  of  isotherms  reveals  a  strong 
convective  motion  in  the  cavity.  The  maximum  heat  transfer,  indicated  by  closely  spaced  isotherms, 
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c.  Ra  =  6350,  R"  =  2.6,y'  =  -2.0;  w,  =  1.90,  by  =  0.238.  d.  Ra„  =  1 7250,  R"  =  2.6,y’  =  -0.77;  w  =  6.70,  V  .  =  -0.75 

n  '  itioa  '  n  ■  iiioa  mm 

and  by  =  1.242. 


Figure  44.  Variation  of  isotherms  and  streamlines  as  function  of  Ra„,  R"  and  y '  (Nguyen  et  al.  1982). 
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occurs  at  the  top  of  the  cavity  for  the  outer  cylinder  and  at  the  bottom  of  the  inner  one.  Figure  44b  has 
the  same  parameters  of  Ran,  R",  but  y '  is  changed  from  0  to  -1  (this  is  equivalent  to  Fig.  38a  and  38b). 
The  flow  pattern  is  a  direct  consequence  of  the  maximum  density  of  water  at  4°C  (the  heavy  dashed 
line  represents  the  4°C  isotherm  and  this  define  the  regions  of  maximum  density).  The  convective 
motion  inside  the  cavity  is  now  considerably  lessened  as  compared  with  the  case  of  y '  =  0. 

Figure  44c  depicts  the  case  of  y'  =  -2  (this  is  the  case  when  T.  =  Tm).  The  local  heat  transfer  is 
maximal  at  the  bottom  of  the  outer  cylinder  and  at  the  top  of  the  inner  one.  This  situation  is 
completely  opposite  the  case  described  in  Figure  44a.  For  the  limited  range  of  -0.75  <  y'  <  =  0.90, 
the  tendency  for  the  clockwise  circulation  pattern  to  form  two  cells  is  shown  in  Figure  44d.  The 
distortion  of  the  isotherm  pattern  in  the  upper  part  of  the  cavity  is  the  result  of  the  intensive 
convection  generated  by  the  clockwise  vortex  located  in  the  region.  This  is  somewhat  similar  to  the 
flow  patterns  exhibited  in  Figure  34c. 

Figure  45a  shows  the  variation  of  dimensionless  temperature  0  as  a  function  of  dimensionless 
radius  R  for  Ran  =  6350,  R"  =  2.6  and  y'  =  0.  The  temperature  distribution  is  very  similar  to  those 
of  the  experimental  values  shown  in  Figure  36a.  Figure  45b  is  the  corresponding  local  Nusselt 
number  distribution  as  a  function  of  angular  coordinate.  The  Nuloc  0  reaches  a  minimum  at  1 80° 
while  the  inner  cylinder  has  a  minimum  heat  transfer  at  0°  (this  is  nearly  identical  to  the  distributions 
shown  in  Fig.  36b). 

The  values  of  Nu .  .  and  Nu, _ are  defined  as 

Joc,i  Joc,o 


M<ioc,  i  =  -  In  R' 


R  + 


39 


3  R 


R+  =  =!/(/?"-  l) 


(64) 


and 


^(oc,o  —  In  R 


R 


+  30 
3/?  + 


R +  =R+  =  R"/(R"-l) 


An  overall  average  A hi  can  be  expressed  as 


,  +  30 
3  R  + 


R+=i/{R"-i) 


In  R' 


,+  30 
3  R  + 


R+=R"/(r"~  l) 


(65) 


(66) 


Figure  46a  shows  the  effect  of  nonJinear_Rayleigh  number  on  Nu  as  a  function  of  y '  for  the  case 
of  R"  =  2.6,  where  Nu  and  the  minimum  Nu  values  increase  as  Ran  increases.  For  all  the  Ran  range 
studied,  the  minimum  Nu  occurred  aty'  values  somewhat  greater  than-1  (=  -0.75).  The  influence 
of  R"  on  Nu  is  shown  in  Figure  46b  for  the  case  Ran  =  1 04.  It  is  clearly  shown  that  Nu  occurs  at  values 
of  y '  >  -1  but  approaches  y '  =  -1  as  R"  increases. 


In  a  square  enclosure 

Watson  (1972)  performed  the  first  analysis  on  the  effect  of  inversion  temperature  on  the 
convection  of  water  in  an  enclosed  square  cavity.  In  this  study,  the  effect  of  temperature-dependent 
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a.  Temperature  profile. 


0* 


ISO* 


b.  Distribution  of  local  Nusselt  number. 

Figure  45.  Variation  of  dimensionless  temperature  9  as  function  of  dimensionless 
radius,  R for  Ra„  =  6350,  R"  =  2.6  and  y'  =  0  (Nguyen  et  al.  1982). 
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b.  For  Ra„  =  1  x  10*  with  R"  as  parameter. 

Figure  46.  Variation  of  overall  Nusselt  number  Nu  with  y'  (Nguyen  et  at.  1982). 


viscosity  was  found  to  result  in  changes  in  magnitude  rather  than  the  characteristics  of  the  flow.  The 
decrease  in  viscosity  in  the  vicinity  of  the  hotter  wall,  accompanied  by  higher  fluid  velocities,  gave 
rise  to  the  increases  in  heat  transfer  indicated  by  the  dashed  curve  in  Figure  47.  It  can  be  seen  that  the 
Boussinesq  model  displayed  a  completely  different  heat  transfer  phenomenon  in  a  system  containing 
a  density  extremum. 

Lin  and  Nansteel  (1987a)  analyzed  the  problem  of  natural  convective  heat  transfer  in  a  square 
enclosure  that  contained  water  near  its  density  maximum.  The  square  enclosure  had  two  sides  of  length 
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Figure  47.  Comparison  of  Nu  from  Boussinesq 
model,  full  equation  and  constant  viscosity  model 
as  function  ofTj  (Watson  1972). 


Figure  48.  Schematics  of  the  square  enclosure  with  thermal 
conditions  (Lin  andNansteel  1987a). 


L  with  the  two  opposing  vertical  walls  maintained  at  Th  and  T,  for  the  temperature  range  0°  <  7,  <  Th 
£  20°C.  The  horizontal  walls  were  insulated  (see  Fig.  48).  The  enclosure  was  supposed  to  be 
sufficiently  deep  to  ensure  an  essentially  identical  flow  field  in  planes  of  different  depth.  In  addition, 
the  following  assumptions  were  made:  1)  the  flow  is  laminar  and  two  dimensional,  2)  physical 
properties,  except  for  the  density  in  the  buoyancy  force  term,  are  constant  and  evaluated  at  the  cold 
wall  temperature  Tc,  and  3)  the  viscous  dissipation  can  be  neglected. 
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The  dimensionless  governing  equations  are 


^  q  |  0-/?|^2  ( Q-R )  +  V2^ 

DlPr  dX 


(67) 


D0  _  J_  v2e 

Dx  Pr 
V2v  = 

The  initial  and  boundary  conditions  are 


(68) 

(69) 


U  =  V  =  0,  0=1,  at  x  =  0 
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e(o,y)  =  i,  e(t,y)  =  o,  at  x  >o 
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^■  =  0  on  the  boundary  with  the  dimensionless  variables  defined  as 

dx 

dr 

x  = 

r  =  y. 

L 

L 

u  =  “L, 

V  =  Vi*- , 

1 

K 

II 

CD 

V 

V 

Th-Tc 

V 

V 
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(70) 


in  which  y  is  the  stream  function,  and  \  is  the  vorticity.  Additional  dimensionless  parameters  ap¬ 
pearing  in  eq  67  and  68  are  the  Rayleigh  and  Prandtl  numbers: 


Ra"  =  ffPm  al  1  (lh  ~  Tc)?  Pr=  X- 
pc  va  ’  a 

in  which  a,  and  q  are  from  the  density-temperature  relation 

P  =  Pm[l-ai|7'-7'ml<?]  (71) 

with  pm  =  999.9720  kg  nr3,  a,  =  9.297 173  x  10^(°CH,  Tm  =  4.029325°C  and  q  =  1 .894816.  The 
density  distribution  parameter  is  defined  as 

R'  =  Is&zJji  (72) 

n-rc 
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which  is  equivalent  to  -1/2  y'defined  in  eq  63.  In  the  case  R'<  0  (i.e.,  for  the  condition  of  Th  >  Tc> 
T  ),  the  fluid  density  increases  monotonically  with  X  across  the  enclosure.  The  density  distribution 
results  in  a  clockwise  circulation  pattern.  When  R'  >  1  (i.e.,  for  Tm>T^>  Tc),  the  distribution  is  re¬ 
versed,  and  hence  a  counterclockwise  pattern  results.  When  R'  is  in  the  range  of  0  <  R'  <  1,  density 
increases  to  p  =  p  at  T  =  T  and  then  decreases.  Hence,  this  heavier  fluid  in  the  enclosure  interior 
descends  while  the  lighter  fluid  adjacent  to  vertical  boundaries  ascends,  giving  rise  to  a  pair  of 
counterrotating  vortices  arranged  horizontally  in  the  enclosure. 

Streamlines  and  isotherms 

For  small  values  of  Ra"\  the  inertial  effects  are  relatively  unimportant,  and  so  an  approximate 
balance  between  viscous  and  buoyancy  forces  results.  Figure  49  shows  a  series  of  dimensionless 
streamlines  for  Ra'"  =  1 03  but  for  a  range  of  R'  values  from  0.4, 0.5, 0.55, 2/3  and  0.75.  At  Ra'"  = 
103,  the  temperature  fields  deviate  from  the  pure  conduction  field;  i.e.,  0  =  1  -X  only  slightly.  However, 
the  convection  effects  on  the  temperature  fields  increase  as  I/?'  -  1/21  increases.  As  R'  increases  (or 
decreases),  the  maximum  density  contour  moves  toward  the  hot  (or  cold)  wall.  As  a  result,  the 
counterclockwise-rotating  cell  on  the  right  for  the  case  of  R'  =  0.4  becomes  stronger  and  larger  at 
the  expense  of  the  clockwise-rotating  cell  on  the  left  as  R'  increases.  For  R'  =  2/3,  the  cell  on  the  left 
was  divided  into  two  separate  clockwise-rotating  cells  in  the  upper  and  lower  left-hand  comers  of 
the  enclosure.  In  this  case,  although  the  maximum  density  contour  (0m = 2/3)  is  locatedapproximately 
at  X=  1  /3,  upward  flow  on  the  hot  walls  occurs  only  near  the  comers  of  the  enclosure.  The  circulation 
of  the  counterclockwise  right-hand  vortex  has  become  so  strong  that  it  drags  (by  virtue  of  the  fluid’s 
viscosity)  relatively  light  fluid  down  along  the  hot  walls.  When  0m  =  R'  =  3/4,  even  the  two  tiny 
comer  cells  are  eliminated  by  the  strong  counterclockwise-rotating  cell  on  the  right,  resulting  in  a 
completely  unicellular  flow.  In  the  instance  R'  =  1,  density  increases  with  temperature  and  hence 
decreases  with  X,  resulting  in  a  single  counterclockwise  rotating  cell. 

Figure  50  shows  the  y  and  0  fields  for  R'  =  1  /2  and  Ra"'  =  1 04  and  1 06.  The  bicellular  flow  structure 
observed  for  Ra"'  =  1 03  and  /?'=  1/2  still  prevails  for  10*  £  Ra"'  £  1 06  with  the  cells  becoming  more 
angular  in  shape.  The  flows  in  the  left  and  right  halves  of  the  enclosure  are  identical  except  for  the 
opposite  direction  of  rotation.  The  large  temperature  gradients  near  X  =  1/2,  Y  £  1  are  due  to  the 
internal  circulation  of  the  two  counterrotating  cells  that  bring  warm  and  cold  water  from  the  hot  and 
cold  walls,  respectively. 

For  a  slight  increase  in  R'  from  R'  =  1/2  to  R'  =  0.55,  there  is  no  spatial  symmetry  (see  Fig.  51) 
because  the  maximum  density  contour  0m = R' = 0.55  has  now  shifted  toward  the  hot  wall.  The  large 
counterclockwise  cell  adjacent  to  the  cold  wall  becomes  more  dominant.  Cold  water  is  swept  across 
the  upper  boundary  and  into  the  upper  left-hand  comer,  which  shifts  the  maximum  density  contour 
0m = 0.55  toward  the  hot  wall  in  the  upper  half  of  the  enclosure.  As  a  result,  there  is  no  upflow  along 
the  upper  portion  of  the  hot  wall  for  Ra'"  =  106  because  the  strong  counterclockwise  vortex  over¬ 
comes  the  (upward)  buoyancy  force  in  the  fluid  directly  adjacent  to  the  upper  portion  of  the  hot  wall . 

Heat  transfer 

The  local  heat  flux  in  the  horizontal  direction  can  be  expressed  as  the  superposition  of  conductive 
and  convective  modes,  i.e., 

q=-k*L+Pocp(T-Tc)u.  (73) 

dx 

Expressed  in  dimensionless  form,  the  Nusselt  number  can  be  rewritten  as 

Nu  ( X,Y )  =  -  qL  .  =  -  ^-  +  PrW .  (74) 

k(Th-Tc)  dX 
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Figure  50.  Dimensionless  stream  function  (y  x  10 J)  and  temperature  (0  =  0(0.1  )J]  contours  for  Ra  //,= 
and  I06  and  R'  =  1/2  (Lin  and  Nansteel  1987a). 

The  average  Nusselt  number  over  a  vertical  cross  section  is 

Nu(X)=  f  Nu(X,Y)dY.  (75) 

Jo 

Figure  52  shows  the  variation  of  A/m  with/?'  for  Ra'"  at  10-\  104, 10<i  and  106.  The  most  striking  feature 
of  this  figure  is  the  minimum  in  heat  transfer  at  R'  =  1/2,  which  is  caused  by  the  symmetric,  dual-cell 
flow  structure  that  results  when  7"h  and  Tc  straddle  the  maximum  density  temperature.  The  dual-cell 
structure  prevents  direct  convective  transfer  between  the  hot  and  cold  walls.  Each  cell  acts  like  an 


Ra'"  =  106 

Figure  51.  Dimensionless  stream  function  (y  x  103)  and  temperature  /0  =  0(0. 1)1]  contours  for  Ra'"=  104 
and  106  and  R'  =  0.55  (Lin  and  Nansteel  1987a). 


insulator,  preventing  warm  (or  cool)  fluid  from  the  hot  (or  cold)  wall  from  coming  in  contact  with 
the  cold  (or  hot)  wall.  The  only  direct  thermal  communication  between  the  two  walls  occurs  at  X  = 
1/2  where  the  warm  and  cold  streams  meet  and  energy  is  transferred  primarily  by  conduction.  Nu 
increases  very  sharply  for  1/2  <R'<  1/2  and  when  Ra'"  is  large,  because  one  of  the  two  cells  in  the 
enclosure  wets  both  walls  (see  Fig.  51). 

Figure  53  shows  the  variation  of  Nu  with  Ra"'  with  /?'  as  a  parameter.  It  can  be  observed  that  as  Ra'" 
becomes  larger,  the  heat  transfer  behavior  is  similar  to  that  observed  for  a  Boussinesq  fluid,  that  is 

Nu  a  (Ra’")0-29 
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independent  of  R The  behavior  for  large  Ra is  seen  to  blend  smoothly  into  the  conduction- 
dominated  regime,  Nu~\,  for  smaller  Ra'". 

In  a  rectangular  enclosure 

Desai  and  Forbes  (1971)  reported  the  first  numerical  investigation  on  free  convection  of  water  (in 
the  vicinity  of  maximum  density)  in  a  rectangular  enclosure.  The  vertical  walls  were  maintained  at 
different  temperatures  and  the  horizontal  walls  insulated.  Two  sets  of  cold  and  hot  wall  temperature 
combinations  (i.e.,  To  =  2°C,  Th  =  6°C  and  Tc  =  0°C  and  7"h  =  8°C)  were  used  for  each  aspect  ratio  of 
1  or  3.  Desai  and  Forbes  (1971)  used  two  representations  of  density-temperature  relations: 


Ps  = 


0.999841+0.000132  sin  fs2l) 

\  8 


and 


pp  =  0.99984247  +  6.460  x  10"5  T-  8.0  x  10-6  T2  . 

Their  finite  difference  approximation  solution  leads  to  the  heat  transfer  across  the  enclosure  as 


in  which 


(where  W  is  the  width  of  the  enclosure),  and  0  is  the  dimensionless  temperature  defined  as 


0"  =  when  T  =(rh  +  rc)/2. 

Th-T 


Table  2  gives  the  steady-state  mean  Nusselt  numbers  as  a  function  of  aspect  ratio  A  (defined  as  height/ 
width),  Te,  Th  and  the  density  temperature  representation.  Clearly,  the  Nu  values  are  higher  for  the 
0°  and  8°C  combination  regardless  of  the  function  of  density  it  is  based  on. 

Nansteel  et  al.  ( 1 987)  theoretically  extended  the  work  of  Lin  and  Nansteel  ( 1 987a)  to  the  case  of 
rectangular  enclosure  at  low  Rayleigh  numbers.  The  schematic  representation  of  the  problem  is  the 
same  as  shown  in  Figure  48  except  for  an  additional  variable,  the  aspect  ratio  A  . 


Table  2.  Steady-state  mean  Nusselt  numbers  (Desai  and 
Forbes  1971). 


1  H/W| 

T  (°C) 

Nu  (based  on  p*) 

Nu  (based  on  pp* 

1 

2 

6 

1.5164 

1.3096 

1 

8 

2.7738 

2.5989 

3 

2 

6 

1.1314 

1.0179 

3 

0 

8 

2.1439 

1.9733 

*  Density  based  on  sine  and  polynomial  function,  respectively. 
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With  the  assumption  of  steady,  laminar  and  two-dimensional  flow  and  constant  physical  properties 
(except  for  the  density  in  the  buoyancy  force  term)  evaluated  at  the  extremum  temperature  and  with 
the  density-temperature  relation  represented  by  eq  2,  i.e., 

P  =  Pm  [l  ~  Y  {T-  T’m)2] 

the  dimensionless  equations  of  motion  and  energy  transfer  are 

w(£  M  -  M  ■ pr  h  -  2*a(9-*'>D  <76> 

and 

dtjt  90^  _  9tjr  90^  _  y2  q  (77) 

BY  dX  dX  BY 

in  which  the  Rayleigh  number  Ra  and  the  density  distribution  parameter  are  given  as 

Ra  =  gyL^Th-Tcf/va 
and 

R  -  (Ym  “  Tc)  /  (Th  —  Tg) . 

The  solutions  of  eq  76  and  77  are  assumed  to  be  expanded  in  powers  of  Ra  and  0,  i.e.. 


V  =  \|/i  Ra  +  \|f2  Ra2  + . . . . 

(78) 

0  =  Rq  +  02  Rq^  +  •  • . . 

(79) 

The  average  heat  flux  at  the  hot  wall  is 


q-H 


H 


~(0»y)  dy 
Bx 


If  expressed  in  terms  of  Nusselt  numbers,  this  becomes 


)’A 

—  (0,  Y)dY. 
o  dX 


Substituting  eq  79  into  eq  80  yields 


Nu  =  Nuq  +  Nu  i  Ra  +  Nu2  Ra 1  +  . 
where 


(80) 


(81) 
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Nut  =  -  -L  *L(0,  Y)dY. 

A  J0  BX 

The  leading  term  NuQ  =  1  is  the  contribution  to  the  heat  transfer  from  the  conduction  temperature 
field  80  =  1  -X.  The  first  correction  Nu]  =  0  is  necessary  because  of  the  antisymmetry  of  0,  with 
respect  to  Y  =  A/2.  Therefore  the  first  non-zero  correction  to  heat  transfer  is 

■A 

nu2  =  -±- 

a  Jo  dx 

which  can  be  determined  in  terms  of  8j  and  x/j  by 

Nu2  =  .  (82) 

Flow  patterns 

Figure  54  shows  the  variation  of  ij/j  with  R'  =  1/2,  for  A  =  1/2,  and  2.  The  strength  of  the 
counterrotating  flow  appears  to  decrease  with  decreasing  A.  Circulation  strength  is  effectively 
quenched  with  decreasing  A  because  of  the  chocking  effect  of  the  horizontal  boundaries  as  they 
approach  one  another.  Figure  55  shows  the  variation  of  \|/,  with  R'  =  0.65  (a  slight  increase  from  Fig. 
54)  and  A-  1/2  and  2.  It  indicates  a  transition  from  a  two-cell  to  three-cell  structure  with  increasing 
A .  This  is  because  the  large  counterclockwise  cell  on  the  right  is  close  to  the  hot  wall  and  able  to  drag 
the  warm  and  lighter  fluid  downward,  overcoming  the  upward  buoyancy  force. 


a.  A  =  112.  b.A  =  2. 

Figure54.  Variation  of  dimensionless  stream  function  (yj  xlO^/forR'  =  l/2(Nansteeletal.  1987). 
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a.  A  =  1 !2. 


b.A  =  2. 


Figure  55.  Variation  of  dimensionless  stream  function  (yjXl&iforR'  =  0.65  (Nansteel  etal.  1987). 


a.  R'  =  112.  b.  R'  =  0.65. 

Figure  56.  Dimensionless  first  correction  0 1  as  function  ofR'  with  A  =  1  ( Nansteel  et 
al.  1987). 


Heat  transfer 

Figure  56  shows  the  first  correction  0j  to  the  conduction  temperature  field  0  =  1  -X  for  a  square 
enclosure.  Figure  56a,  for/?'  =  1/2,  clearly  shows  that  0j  is  antisymmetric  with  respect  to  the  en¬ 
closure  midplanes  X=\/2  and  Y = All.  As  a  consequence,  9{  is  a  positive  (or  negative)  perturbation 
on  0Q  in  the  upper  left-hand  (or  right  -hand)  and  lower  right-hand  ( or  left-hand)  quadrants  of  the 
enclosure.  Hence,  isotherms  (0)  will  be  slightly  compressed  toward  the  vertical  midplane,  X  =  1/2, 
in  the  upper  half  of  the  enclosure  and  displaced  outward  from  the  midplane  in  the  lower  half.  Figure 
56b  shows  0 1  contours  for  the  case  of/?'  =  0.65.  It  indicates  that  0 ,  results  in  a  negative  perturbation 
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Figure  57.  First  correction  to  heat 
transfer,  Nu,  as  function  of  R' 
(Nansteel  etal.  1987). 


R 


to  the  temperature  field  in  the  upper  half  of  the  enclosure  and  a  positive  perturbation  in  the  lower 
half.  This  results  in  isotherms  being  displaced  slightly  toward  the  hot  wall  in  the  enclosure’s  upper 
half  and  toward  the  cold  wall  in  the  lower  half.  Figure  57  shows  the  variation  of  the  first  correction 
to  the  heat  transfer  Nu2  with  R'  for  A  =  1/2, 1 ,  and  2.  Nu2  reaches  a  minimum  at  R'  =  1/2.  It  should 
be  noted  that  the  perturbation,  Nu2,  to  the  pure  conduction  heat  transfer  rate  increases  rapidly  and 
symmetrically  with  respect  to  R'  =  1/2.  The  magnitude  of  Nu2  as  well  as  the  sharpness  of  the 
minimum  at /?'=  1/2  is  clearly  seen  to  be  an  increasing  function  of  A.  For  fixed/?'  ^  1/2  with  increasing 
A,  there  may  be  radical  change  in  the  basic  structure  of  the  flow.  As  indicated  in  Figures  54a  and  55b, 
the  flow  changes  from  a  dual-cell  to  tri-cell  structure  in  which  the  large  counterclockwise  right-hand 
cell  wets  both  the  hot  and  cold  walls.  This  direct  contact  of  cells  with  both  vertical  boundaries  greatly 
enhances  the  convective  contribution  to  the  cross-cavity  energy  transfer. 

In  a  horizontal  layer 

Forbes  and  Cooper  (1975)  presented  the  first  analytical  study  on  natural  convection,  caused  by 
cooling  from  above,  of  a  horizontal  water  layer  initially  at  temperatures  of  either  4°C  or  8°C.  The 
water  was  confined  laterally  and  underneath  by  rigid  insulators  and  the  upper  horizontal  surface  was 
subjected  to  1 )  a  constant  0°C  temperature  and  a  rigid  conducting  boundary,  and  2)  a  free  water-to- 
air  convection  boundary  condition,  in  which  the  convective  heat  transfer  was  held  constant  at  values 
of  5.68  or  284  W/m2K  while  the  temperature  of  the  ambient  air  is  kept  at  0°C.  In  their  study,  the  ratios 
of  width  to  depth  (W/D)  were  1, 3  and  6  (W/D  is  the  reciprocal  of  aspect  ratio  A). 

If  the  flow  is  considered  to  be  laminar  and  two-dimensional  (this  restricts  the  dimension  of  D  to 
the  order  of  1  cm),  and  density  variations  are  considered  only  in  the  buoyancy  term,  the  governing 
dimensionless  equations  of  motion  and  energy  can  be  written  as 

Vorticity 


dSL+LL*L  +  \L*L  =  ±(Ra*)+  V2£2 
at  Pr  dx  Pr  dYdX 


(83) 


and 
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Energy 


Pr  — +  U  —  +  V  — =  v2e 
5t  ax  ar 

Stream  function 

Q  =  -  V2  v 

Velocities 

U  =  v/  =  _ 

ay  ’  ax ' 


(84) 


(85) 


(86) 


The  above  equations  have  been  made  dimensionless  through  the  use  of  the  following  dimensionless 
variables: 


(/  =  ,  V  =  v&  ,  t  =  VI  (87) 

a  a  q2 

¥-*.  o-El. 

a  a 

The  initial  and  boundary  conditions  may  be  written  as 

iir(x,y,o)  =  o,  e(x,y,o)=i 
f/  (x,y,o)  =  v  (x,y,o)=o 

U  =  V  =  0  on  all  boundaries 
and 

—  (o,y,  t)  =  y,t )  =  —(x,  o, t)  =  o. 
ax  dY\D  I  BY 

At  the  free  surface  (Y  =  1)  the  boundary  condition  is  either: 

0""  (X,  1,  x)  =  0  where  0""  =  T~Ts  (88) 

T0-Ts 


for  fixed  surface  temperature  or 


If  0'"  is  defined  as 


(90) 


Ra*  =  Prp[p-p(ew)]gP3 

H2 

where  p  is  the  density  at  (Ts+Tq)P.  in  which  Tq  and  T$  are  the  initial  water  and  water  surface  tem¬ 
perature,  respectively.  The  temperature-density  relation  is  expressed  by 

p  (T)  =  0.99984247  +  6.46  x  10-5  T  -  8.0  x  10-6  T2. 


Figure  58.  Transient  dimensionless  stream  lines:  W/D  =  1.0.  T0  =  8°C,  t  =  055  s 
(x  =  0.002),  T S~0°C  (Forbes  and  Cooper  1975). 
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a.  Transient  streamlines. 


Figure  59.  Transient  dimensionless  streamlines  and  isotherms:  W/D  =  1.0,To=8°C, 
t  =  250  s  (x=  1.0),  and  Ts  =  0°C  (Forbes  and  Cooper  1975). 


Streamlines  and  isotherms 

The  implicit  alternating  direction  method  was  used  to  solve  the  governing  equations,  and  the 
results  are  expressed  in  transient  streamlines  and  isotherms.  Figures  58-6 1  show  a  series  of  transient 
streamlines  and  the  corresponding  isotherms  for  W/D  =  1 .0,  T0  =  8°C,  and  a  constant  surface  tem¬ 
perature  of  Ts = 0°C,  with  time  values  of  0.5  s  (x = 0.002),  250  s,  and  62.5  min.  The  streamlines  shown 
in  Figure  58  are  identical  to  the  flow  configurations  observed  in  both  experimental  and  theoretical 
investigations  of  natural  convection  in  horizontal  layers  of  fluid  heated  from  below.  In  this  figure, 
the  4°C  isotherm  lies  very  near  the  water  surface,  and  the  entire  cavity  region  is  hydrodynamically 
unstable.  The  water  current  ascends  in  the  central  region  between  the  two  eddies,  and  descends 
adjacent  to  the  two  lateral  sides  of  the  cavity. 

At  x  =  0. 1  or  250  s,  (Fig.  59a  and  b),  the  eddy  that  appeared  on  the  left  side  of  the  flow  field  at  x 
=  0.002  has  momentarily  disintegrated.  The  heavy  dotted  line  shows  the  4°C  isotherm;  this  region 
of  maximum  density  provides  the  driving  force  for  the  eddy  motion.  As  the  region  of  the  maximum 
density  descends  further,  the  4°C  isotherm  separates  the  two  eddies  as  shown  in  Figure  60a.  The 
eddy  entirely  beneath  the  4°C  isotherm  is  therefore  the  result  of  unstable  density  gradients. 

This  circulation  flow  exerts  a  viscous  shearing  stress  upon  the  stable  water  above  the  region  of 
maximum  density,  and  this  shearing  stress  provides  the  energy  necessary  for  eddy  formation  in  the 
region  of  hydrodynamic  stability.  The  straight  parallel  isotherms  that  appear  in  the  upper  portion  of 
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- 7 


b.  Transient  isotherms. 
Figure  59  (cont’d). 


Table  3a.  Maximum  and  minimum  values  of 
stream  functions,  with  a  constant  temperature 
boundary  (Desai  and  Forbes  1971). 


Figure 

W/D 

T 

o 

w 

T  max 

Ay 

58 

1 

8 

3.0 

-3.0 

1.0 

59(a) 

1 

8 

0.0 

-60 

0.40 

60(a) 

1 

8 

0.5 

-0.4 

0.05 

61(a) 

1 

8 

O.OOI28 

0.0 

0.00008 

Table  3b.  Maximum  and  minimum  dimension* 
less  temperatures,  with  a  constant  temperature 

boundary  (Desai  and  Forbes  1971). 

Figure  W/D  T. 

n 

//// 

®max 

tw 

e  . 

mm 

AO"" 

59(b)  1  8 

0.72 

0.0 

0.04 

60(b)  1  8 

0.60 

0.0 

0.03 

61(b)  1  8 

0.034 

0.0 

0.002 
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a.  Transient  streamlines. 

F inure  60.  Transient  dimensionless  streamlines  and  isotherms:  W/D  =  1.0,  To  =  H  °C, 
t  =  500  s  (t  =  2.0)  and  Tf  =  0°C  (Forhes  and  Cooper  1975). 

Figure  60b  indicates  that  the  heat  transfer  in  the  stable  region  is  largely  by  conduction.  The  eddy  in 
the  unstable  region  continues  todistort  the  temperature  distribution  and,  therefore,  heat  is  transferred 
in  this  region  primarily  by  convection.  Figures  6 1  a  and  b  depict  the  flow  pattern  and  the  temperature 
distribution,  respectively,  at  T  =  5.0  (or  62.5  min).  The  4°C  region  has  now  passed  out  of  the  domain 
of  the  system  and,  therefore,  the  fluid  body  is  entirely  stable.  The  large  eddy  appearing  in  Figure  6 1  a 
is  the  same  eddy  that  previously  appeared  above  the  region  of  maximum  density  (Fig.  60a).  The 
momentum  of  this  eddy,  at  this  time,  has  not  yet  been  dissipated  by  the  opposing  viscous  forces. 
Consideration  of  the  magnitudes  of  the  stream  function  (see  Table  3a,  3b)  reveals  that  the  fluid 
velocities  are  much  smaller  than  those  depicted  in  Figure  60a.  Indeed,  the  fluid  motion  is  so  minute 
that  the  temperature  field  appears  undisturbed  (Fig.  6 1  b)  and  the  water  seems  to  have  stratified.  Heat 
transfer,  therefore,  at  these  times  is  purely  conductive. 

In  the  case  of  Ta  =  4°C,  the  imposition  of  0°C  upper  surface  boundary  temperature  cannot  lead  to 
eddy  formation,  as  the  density  gradient  in  this  case  is  always  stable.  The  fluid  motion  is  not  sufficient 
to  cause  distortion  of  the  temperature  field  and  the  heat  transfer  is  purely  conductive. 

For  the  case  of  Tq  =  8°C,  numerical  results  obtained  for  the  W/D  =  3  and  6  enclosures  were  es¬ 
sentially  identical  to  those  for  W/D  =  1,  differing  only  in  the  number  of  eddies  produced  by  the 
motion.  Three  and  seven  pairs  of  coupled  eddies  were  produced  in  W/D  =  3  and  6,  respectively,  and 
the  4°C  isotherm  descended  through  the  layer.  In  all  cases  studied,  the  motion  was  "damped  out" 
as  the  4°C  isotherm  reached  the  bottom,  with  a  stable  water  layer  resulting  for  the  entire  cavity. 


b.  Transient  isotherms. 


Figure  60  (cont’d). 


Convective  boundary 

For  the  case  of  h  =  5.68  W/m2K,  Forbes  and  Cooper  (1975)  reported  that  it  took  much  longer  to 
set  up  the  region  of  maximum  density.  The  following  results  will  summarize  the  work  on  h  =284  W/ 
m2K.  As  indicated  eq  90,  i.e.. 


dQ'" 

dY 


(0'")y=i  - 


The  value  of  (0"')y=|  =  1.0  for  both  values  of  h.  Since  D  and  k  are  constants,  a  50-fold  increase  in  the 
magnitude  in  h  will  produce  50-fold  increase  in  the  magnitude  of  the  initial  temperature  gradient  at 
the  water  surface.  Figures  62a  and  b,  63a  and  b,  64a  and  b  show  the  streamline  and  isotherms  for  W/ 
D  =  3,Tq  =  8°C,  h  =  284  W/m2K  for  a  real  time  of  3, 4  and  5  minutes,  respectively.  Figures  65a  and 
b  show  the  case  of  W/D  =  6.0,  7q  =  8°C,  h  =  284  W/m2K  with  t  =  4  min.  Tables  4a  and  b  show  the 
maximum  and  minimum  values  of  stream  functions  and  dimensionless  temperature,  respectively. 
Figures  62a  and  b  reveal  the  flow  pattern  and  temperature  distribution  3  min  after  the  initiation.  Since 
the  region  of  maximum  density  is  very  near  the  upper  surface  at  this  time,  the  flow  pattern  shown  in 
Figure  62a  is  very  similar  to  the  flow  that  occurs  in  horizontal  fluid  layers  heated  from  below.  At  t  = 
4  minutes  after  initiation,  the  4°C  isotherm  descends  far  enough  through  the  liquid  to  produce  two 
layers  of  eddies.  The  lower  eddies  are  the  result  of  an  unstable  density  gradient  and  the  upper  eddies 
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Figure  61.  Transient  dimensionless  streamlinesandisotherms:  W/D  =  l  .0,  Tc  =  <5CC, 
t  =  62.5  min  (x  =  5.0)  and  T S  =  0°C  ( Forbes  and  Cooper  1975). 


Table  4a.  Maximum  and  minimum  values  of 
stream  functions,  with  a  constant  temperature 
boundary  (Desai  and  Forbes  1971). 


Figure 
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-2.4 

63(a) 

3 

-0.22 

0.02 

64(a) 

3 

0.01 

65(a) 

6 

-0.22 

0.02 

Table  4b.  Maximum  and  minimum 
dimensionless  temperatures,  with  a 
constant  temperature  boundary 
(Desai  and  Forbes  1971). 
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b.  Transient  isotherms. 
Figure  61  (cont’d). 


a.  Transient  streamlines. 

Figure  62.  Transient  dimensionless  streamlines  and  isotherms:  W/D  =  3.0,  T 0  =  8°C, 
t  =  min,  and  h  =  284  Wln^k  ( Forbes  and  Cooper  1975). 
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b.  Transient  isotherms. 


Figure  62  (cont'd).  Transient  dimensionless  streamlines  and  isotherms:  W/D  =  3.0,  T0 
=  8°C,  t  =  min,  and  h  =  284  W/rt?k  (Forbes  and  Cooper  1975). 


a.  Transient  streamlines. 


Figure  63.  Transient  dimensionless  streamlines  and  isotherms:  W/D  =3.0,To=8  °C, 
t  =  4  min,  and  h  =  284  Wlm?K  (Forbes  and  Cooper  1975). 


are  caused  by  the  action  of  viscous  shearing  stress  upon  the  water  above  the  unstable  region.  Figures 
64a  and  b  provide  an  even  clearer  picture  of  the  effect  of  the  region  of  maximum  density  on  free 
convection  in  water.  The  4°C  isothermal  line  is  observed  to  exist  near  the  bottom  of  the  enclosure 
and  the  dimensionless  temperature  gradient  ( dB'"/dY)  is  almost  constant  in  the  hydrodynamically 
stable  portion  of  the  water.  Figures  65a  and  b  show  the  case  for  W/D  =  6.  The  transient  formation 
of  streamlines  and  dimensionless  temperature  distribution  are  very  similar  in  nature  as  in  the  case 
of  W/D  =  3.  However,  it  contains  seven  pairs  of  eddies  instead.  No  conclusive  reason  is  given  for 
the  asymmetrical  nature  of  the  streamline  and  the  temperature  distribution. 
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a.  Transient  streamlines. 


Figure  64  .  Transient  dimensionless  streamlines  and  isotherms:  W/D=  3.0,  Te  =  8°C, 
T  =  5  min,  and  h  =  284  Wln^K  (Forbes  and  Cooper  1975). 


Figure  65.  Transient  dimensionless  streamlines  and  isotherms:  W/D  =  6.0,  T0  =  8°C, 
t  =  4  min,  and  h  =  284  WlmrK  (Forbes  and  Cooper  1975). 


In  a  circular  and  confined  melt  layer 

Yen  (1968)  and  Yen  and  Galea  (1969)  were  the  first  to  investigate  the  natural  convection  and  the 
temperature  distribution  caused  by  the  density  inversion  of  the  continuously  deepening  layer 
(originally  there  is  only  a  single  phase,  i.e.,  a  solid  ice  melt  layer  forming  upon  the  application  of  a 
constant  surface  temperature).  In  their  study,  both  modes  of  melting  were  employed  (i.e.,  melting  from 
below  and  the  top).  They  found  that  the  onset  of  convection  and  the  transition  from  a  conductive  to 
convecti  ve  heat  transfer  mode  depends  on  the  imposed  boundary  temperatures,  a  rather  strong  function 


71 


ICE 


b)  Melting  from  Above 


Figure  66.  Schematic  representation  of  the  interdependence  of 
stable  and  unstable  regions  with  the  thermal  boundary  conditions. 


of  temperature  in  melting  from  below  but  a  considerably  weaker  one  in  the  case  of  melting  from 
above.  Figure  66  is  a  schematic  representation  of  the  interdependence  of  stable  and  unstable  regions 
with  the  thermal  boundary  conditions.  Figures  67  and  68  show  some  typical  temperature  distributions 
in  the  melt  layer  at  various  times  after  initiation  of  the  melting  (Yen  1 984).  It  can  be  seen  that  the 
time  elapsed  to  deviate  from  linear  distribution  depends  on  the  boundary  temperature  imposed.  For 
melting  from  below,  the  higher  the  boundary  temperature  the  shorter  the  time  needed  to  change  from 
conduction  to  convection.  On  the  other  hand,  for  melting  from  above,  since  the  maximum  density 
region  is  near  the  water/ice  interface,  the  overlying  stable  layer  becomes  thicker  as  upper  boundary 
temperatures  increase,  extending  the  4°C  region  upward  and  thus  reducing  its  effectiveness  in 
overturning  the  water  near  the  ice  surface.  The  major  difference  in  the  temperature  distribution  in 
the  pseudo-steady  state  is  that  the  constant  temperature  zone  depends  on  the  boundary  temperature 
in  melting  from  below,  while  the  temperature  in  the  constant  temperature  region  is  ~3.2°C  inde¬ 
pendent  of  the  upper  boundary  temperature  in  melting  from  above  (see  Fig.  68c). 

Yen  (1980)  derived  heat  flux  expressions  for  the  melt  layer  formed  by  melting  from  both  below 
and  above.  In  his  experiments  where  the  melting  rate  was  determined,  the  total  upward  or  downward 
heat  flux  for  melting  from  below  or  above  was  evaluated  by 

q=  p  &  [Lf+  C,  (Tm -  Ti0)]  +  ^  [z  (r) (T  +  Tm)]  ~kM\  (91) 

dt  2t  dz  lz  =  n. 
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Mean  Water  Temperature  (°C) 


Water  Layer  Depth  (cm) 


Water  Layer  Depth  (cm) 


a.  Before  onset  of  convection. 


b.  Development  of  the  convective  layer. 


c.  Pseudo-steady  state. 


Figure  67.  Mean  temperature  profiles  for  melting  from  below  (after  Yen  1984). 


73 


Mean  Water  Temperature  (°C) 


T2("C)  T0  (  C ! 
to)  11.92  -11.5 

(a)  1798  -11.5 

(•)  27.23  -11.5 

(A)  32.14  -11.5 

(*)  39.75  -11.5 


vt  =  206  (min.) 


T2(°C)  T0(°C) 
(o)  1 1.92  -|  1.5 

(o)  17.98  -11.5 

(•>  27.23  -II. 
(A)  32.14  -II 
(»)  39.75  -II. 


-  238  (min.) 


I  2  3 

Water  Layer  Depth  (cm) 


2  3 

Water  Layer  Depth  (cm) 


a.  Before  onset  of  convection. 


b.  Development  of  the  convective  layer. 


\t  =  651  (min.) 


T2(°C)  Tq  ( 0C ) 

( o )  n.92  -11.5 

(□)  17.98  *115 

(a)  39.75  -11.5 


2  4  ~6 

Water  Layer  Depth  (cm) 

c.  Pseudo-steady  state. 


Figure  68.  Mean  temperature  profiles  for  melting  from  above  (after  Yen  1984). 
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where  L{  =  latent  heat  of  fusion 
C j  and  Tj  =  specific  heat  and  temperature  of  ice 
p  =  water  density 
rj0  =  the  ice  initial  temperature. 

The  term  (pc  J2t)[z(t)(Tm+Tj\  represents  the  mean  sensible  heat  content  variation  of  the  entire  layer 
of  a  depth  z(t).  Since  Tm  =  0°C,  therefore,  T  stands  for  either  Tl  or  T2.  The  contribution  of  this  term 
to  the  overall  heat  flux  was  found  to  be  much  more  significant  in  the  case  of  melting  from  above 
(since  the  melting  rate  is  much  slower  in  the  case).  For  experiments  in  which  the  mean  layer 
temperature  was  measured  as  the  melting  progressed,  the  heat  flux  was  approximated  by 


(92) 


in  which  ( dT /dz)2  and  (dT /dz) ,  are  the  mean  temperature  gradients  at  the  stable  region  near  the  upper 
boundary  for  melting  from  above,  and  at  the  lower  warm  plate  for  melting  from  below.  Subscripts 
1  and  2  indicate  the  beginning  and  end  of  each  period  and  At  is  the  total  time  period  for  calculating 
a  value  of  q.  At  least  a  dozen  or  more  of  the  periods  (with  varying  durations)  were  used.  The  heat 
from  above  was  found  to  be  a  rather  weak  function  of  T2  and  can  be  represented  by 

q  =  177  (T2)0303  .  (93) 

Eq  93  is  valid  for  T2  ranging  from  1 1 .75  to  39.9°C.  For  melting  from  below,  the  heat  flux  is  found 
to  be  a  linear  function  of  T{  and  can  be  expressed  as 

q  =-  1900  +  315  (T,)  (94) 

for  7,  ranging  from  7.7  to  25.5°C.  Higher  values  of  Tx  and  thus  higher  temperatures  in  the  convective 
zone,  resulted  in  greater  convective  motion  in  the  unstable  region  and  consequently  reduced  the 
thickness  of  the  stable  layer  adjacent  to  both  the  ice  and  the  thermal  boundary.  The  work  of  Townsend 
( 1 964)  and  Adrian  ( 1 975)  is  similar  to  the  case  of  melting  from  above.  However,  in  their  investigations, 
a  rather  deep  invariant  water  layer  was  used  throughout  the  experiment  and  there  was  no  phase 
transition  taking  place  at  the  bottom  of  the  tank.  Both  Townsend  (1964)  and  Adrian  (1975)  reported 
a  nearly  constant  heat  flux  of  approximately  340  W/m2  independent  of  the  initial  water  temperature. 

Melting  in  bulk  water 

Merk  (1954)  reported  the  most  comprehensive  and  earliest  theoretical  study  of  melting  free 
convective  heat  transfer,  such  as  that  for  spheres  and  cylinders.  Employing  a  third-order  density- 
temperature  polynomial  of  water  density  and  applying  the  Von  Karman-Pohlhausen  integral 
method  for  the  case  of  Pr  »  1 ,  he  successfully  solved  the  boundary  layer  equation  and  developed 
a  general  Nusselt  number  ratio: 


Nil.  = 
Nu0 


P  l 


m^p2m  +  -841Lp3„)  li±MP 
13671  31031  I  1  -St  . 


(95) 


where  S  is  the  shape  factor  of  the  temperature  profile  and  is  connected  to  the  Stefan  number  St  = 

cpVnTTJL  fas 
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or 

St  =  6S  — ,  (97) 

(S+2)2 

where  is  the  bulk  water  temperature.  The  values  of  P^P2  and  P2  are  expressed  as 

P,  =  1  -  J2_s  +  J5_  52 
217  434 

/>2  =  1  —  -265.5+  JQ3_  52__l2_  S3 
1343  2686  2686 

and 

p,  =  1  -  Mils  +  -2522,  $2 _J9£_  53  +  — 1£_  S4 
8471  16942  8471  16942 

The  values  of  m  and  n  are  defined,  respectively,  as 
m  =  (p2  +  3  P3  rJe^/vlL 

«  =  £  e2,/^) 

where©  -T-T  and  N  and  P  are  defined  as 

m  m«»  ^00 

N  =  1  +  P','  T00+  &  Tj+  &  Tl 

and 

V~  =  (pV  +  &  T„+3%rl)/N 

in  which  the  P’s  are  the  coefficient  in  the  formula  for  the  specific  volume  of  the  water,  i.e., 

1  =  _L  (1  +  p "t+%T2+%  T-3)  (97a) 

P  Po 

Nu0  is  the  Nusselt  number  for  the  case  of  St  =  m  =  n  =  0.  Merk  ( 1 954)  reported  that  for  large  values 
of  the  Prandtl  number,  neither  the  shape  of  the  body  nor  die  position  on  the  surface  influences  the 
ratio  Nu/Nu0.  Therefore,  this  ratio  can  be  replaced  by  Nu/Nuq,  which  can  be  expressed  by  the 
following  if  only  the  effect  of  melting  is  considered  (i.e.,  for  m  =  n  =  0) 

_  'i  _  .82.  j  +  JJLs2  ,  1,/4 

Mil.  =  _ 212 - 424 —  (1  +  1  s)5  .  (98) 

Nu0  1-ij  +  lS2  '  2  1 

L  2  4  J 
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Figure  69.  The  influence  of  melting  (S  <  0)  and  solidi¬ 
fication  (S  >0)on  the  heat  transfer  in  thermal  convec¬ 
tion  for  large  Prandtl  numbers  (after  Merk  1954). 


Figure  69  is  a  graphic  representation  of  eq  98;  it  clearly  shows  that  for  S  <  0  or  for  melting,  <  Nuq. 
For  solidification,  i.e.,  S>0,Nu  >  Nuq  For  the  case  of  no  melting,  but  with  convective  inversion, 
then  S  =  0  and  hence  P{  =  P2  =  P3  =  1 .  Eq  95  becomes 

=  (1  +  0.4912  m  +0.2730  n) 1/4  .  (99) 

Nuo 

By  defining  Nuq  for  large  values  of  the  Prandtl  number  as 
_  t,3  \l/4 

Nuq  =  C{GrPi) =  c(®L-  Mm)  000) 

eq  99  can  be  written  as 


f  =  - (&L__  =  [Mm(l  +0.4912  m  +  0.2730  n)]1/4  (101) 

(C^r) 

where /is  a  dimensionless  number.  For  small  values  of  T^,  p^  =  P,"  +  2P2"  7^,  m  =  P2'0m/Poo.  and 
n  *  0,  eq  101  can  be  transformed  to 

f  =  [1.509P2  \T~-  riv|  \Tm-  r-|  ] 1/4  002) 

in  which  T.y  is  the  inversion  temperature: 


7iv  =  -0.663  ^1-  -  0.326  Tm. 

& 


003) 
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Figure  70.  Convective  inversion  for  melting  ice  in 
water  of  temperature  T^.  The  arrows  indicate  the 
direction  of the, flow  along  the  surface.  The  dashed  curve 
shows  the  behavior  neglecting  melting,  while  the  full 
curve  of  the  melting  is  taken  into  account  (after  Merk 
1954). 


Figure  71 .  Comparison  of  theoretical  ( eq  109)  and 
experimental  Nu  for  unidirectional  convection  (after 
Schechter  and  Ishin  1958). 


Figure  7L  Comparison  of  theoretical  (eq  110)  and  experi¬ 
mental  Nu  for  inverted  regime  (after  Schechter  and  Ishin 
1958). 

Using  appropriate  values  of  {$("and  p2"and  taking  Tm = 0°C,  Merk  reported  a  value  of  Tjv  =  5.005°C. 
indicating  that  the  inversion  temperature  for  melting  ice  in  water  is  somewhat  greater  than  =•  4°C. 
The  significance  of  the  inversion  temperature  is  clearly  seen  from  eq  102  at  7^  =  Tjv,  Nu  =  0,  and 
the  direction  of  the  flow  in  the  boundary  layer  along  the  surface  of  the  body  is  inverted  (for  <  Tjv 
the  flow  is  upward  and  for  T 'm  >  Tjv  it  is  downward).  Using  the  general  eq  95,  Merk  derived  the 
minimum  Nusselt  number  at  Tjv  =  5.30°C  for  S  =  0  (no  melting)  and  Tjv  =  5.3 1  °C  for  S  <  0  (with 
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melting).  Figure  70  shows  the  effect  on  the  value  of/(i.e.,  the  value  of  Nu)  with  and  without  melting 
with  the  convective  inversion.  It  is  clearly  shown  that  the  effect  of  melting  is  only  appreciable  for 
T„  >  T.  and  may  be  neglected  for  <  Tiv.  Experimental  results  of  Dumor6  et  al.  (1953)  and 
analytical  results  from  the  non-melting  vertical  plate  study  by  Ede  (1955)  generally  confirmed 
Merk’s  findings. 

Tkachev  ( 1 953),  using  photographic  techniques,  reported  a  minimum  Nusselt  numberformelting 
ice  cylinders  at  5.5°C  and  was  the  first  to  notice  the  peculiar  nature  of  the  maximum  density 
boundary  layer.  He  suggested  that  under  certain  conditions  the  boundary  layer  might  be  split,  with 
a  region  of  predominantly  upward  motion  immediately  adjacent  to  the  ice  surface  and  a  region  of 
downward  motion  outside  this.  Tkachev  conducted  melting  experiments  on  spheres  as  well  as  on 
vertical  and  horizontal  cylinders.  Using  the  same  initial  cylinder  diameters  but  with  various  bulk 
water  temperatures,  he  found  that  the  coefficient  of  heat  transfer  is  lowest  for  a  water  temperature 
of  about  5.5°C.  He  correlated  his  data  with  the  following  dimensionless  expressions  as 

Numd  =  0.40  ( GrPr)l&  (104) 

and 

Numd  =  0.104  (GrPr)md  (105) 

for  cylinders  for  the  values  in  the  range  of  1 02  <  GrPr  <  1 07  and  (GrPr)md  >  >  1 07,  respectively.  The 
corresponding  equations  for  spheres  are 

Numd  =  0.54  (GrPr)\&  (106) 

for  laminar  flow  [103  <  (GrPr)md<  107]  and 

Numd  =  0.135  {GrPr)\&  (107) 

for  turbulent  motion  [(GrPr)md  >  107].  Subscript  md  represents  the  physical  properties,  and  the 
diameter  of  the  sphere  or  cylinder  was  evaluated  at  the  arithmetic  mean  temperature  [i.e.,  (T^T^)/ 
2]  and  mean  diameter  [i.e.,  (dQ+df)/ 2]  where  dQ  and  df  are  the  initial  and  final  diameter  at  the  end  of 
the  experiment.  Based  on  his  experimental  data,  Tkachev  further  presented  an  expression  for 
determining  the  time  required  for  completing  melting  as 

StFomdNumd  =  0.305.  (108) 

The  suggestion  of  split  boundary  layer  flow  was  verified  by  the  analytical  and  experimental  work 
of  Schechter  and  Isbin  ( 1 958),  in  which  an  isothermal,  vertical,  non-melting  plate  was  used.  Figures 
71  and  72  compare  the  theoretical  and  experimental  results  for  the  unidirectional  and  inverted 
convections,  respectively.  The  theoretical  curves  in  Figures  71  and  72  are  given  by 
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Figure  73.  Nusselt  number  as 
function  of  Rayleigh  number  (GrPr) 
(after  Vanier  and  Tien  1970). 


where  m  and  n  are  the  same  as  appeared  in  eq  95.  They  reported  that  a  test  of  the  type  of  region  that 
will  prevail  for  given  conditions  of  plate  and  bulk  temperature  can  be  stated;  i.e.,  if 

.1 .+  m+!L  2  0. 

3  5  7 

there  will  be  normal  convection  (unidirectional),  and  if 

1+  2L+  a.  <  o, 

3  5  7 

there  will  be  inverted  convection.  Schechter  and  Isbin  (1958)  concluded  that  the  heat  transfer 
coefficient  can  be  predicted  for  both  regions  with  a  deviation  in  Nusselt  number  of  ±  1 0%,  provided 
the  absolute  value  of  ( 1/3  +  m/5  +  n/7)  >  0.05  by  use  of  eq  1 09  or  1 1 0,  depending  on  the  convective 
region.  However,  the  boundary  layerequations  as  approached  either  by  the  Von  Karman-Pohlhausen 
integral  method  or  by  the  similarity  transformation  method  did  not  yield  meaningful  results  under 
split-flow  conditions. 

To  resolve  this  problem  Vanier  and  Tien  ( 1 968)  used  an  accurate  numerical  method  to  solve  the 
similarity  equations  for  a  semi-infinite  vertical  plate  at  constant  temperature  7"w  immersed  in  an 
indefinitely  large  volume  of  water  at  bulk  temperature  T oo.  They  reported  that  a  new  solution  is 
necessary  for  every  combination  of  7"w  and  Tm,  By  obtaining  several  hundred  such  solutions,  the 
authors  were  able  to  map  out  temperature  zones  for  each  flow  regime.  The  split  boundary  layer  was 
found  to  be  confined  to  two  distinct  triangular  regions  within  which  the  similarity  equations  become 
quite  intractable.  They  confirm  the  findings  of  Merk  ( 1 954)  that  the  melting  heat  transfer  rates  were 
very  similar  to  those  for  the  case  of  non-melting. 

Vanier  and  Tien  (1970)  conducted  experimental  work  aimed  at  relating  their  numerical  plate 
results  to  a  more  practical  geometry  of  the  sphere  (including  the  effect  of  changing  body 
configuration).  This  was  also  partially  motivated  by  lack  of  detailed  analysis  and  correlations  of  the 
experimental  results  on  the  melting  of  ice  spheres  and  cylinders  presented  by  Dumord  et  al.  (1953) 
and  Tkachev  ( 1953).  They  presented  their  results  in  a  least-squares-fitted  semi-empirical  equation  as 

Nu  =  2+  C(GrPr)m  (111) 

and  found  that  for  Tm  >  7°C,  that  the  results  did  not  appear  to  be  affected  by  the  maximum  density  and 
that  the  most  appropriate  value  of  C  is  0.422  ±  0.006  in  the  range  of  1.7xl06<  GrPr  <  2.4x19®  (Fig. 
73,  curve  b).  However,  for  Tm  <  7°C,  nearly  the  same  value  of  C  is  found  but  with  considerably  more 
scattering  (twice  the  standard  deviation),  indicating  the  need  at  least  to  incorporate  another  parameter 
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Figure  74.  Nusselt  number  variation 
with  bulk  temperature  for  approxi¬ 
mately  constant  sphere  diameters 
( after  Vanier  and  Tien  1970). 


Figure  75.  Correspondence  between 
the  melting  of  flat  plates  and  spheres 
of  ice  ( after  Vanier  and  Tien  1970). 


such  as  T^  to  adequately  describe  the  heat  transfer  under  these  conditions.  This  was  done  by  separating 
the  low  temperature  data  into  positive  and  negative  deviations.  For  Tm  <  3.8°C,  Nu  values  are  higher 
than  expected  (curve  a),  while  for  4. 1°  <  <  7. PC,  Nu  values  are  too  low  (curve  c).  To  check  the 
one-quarter  power  assumption  in  eq  1 1 1 ,  a  two-parameter  fit  was  carried  out,  resulting  in 

Nu  =  2  +  0.437  (GrPr) 0-248  (112) 

which  provides  an  excellent  verification.  However,  when  the  Grashof  number  was  calculated  by  using 
an  arithmetic  mean  temperature  basis  of  Tm,  the  constant  C  was  found  to  be  0.52  for  T^  >  1 4°C.  This 
provides  a  remarkable  agreement  with  the  results  reported  by  Tkachev  (1953).  The  effect  of  sphere 
diameter  and  maximum  density  on  heat  transfer  can  be  seen  in  Figure  74  where  the  curves  are  in  general 
agreement  with  the  flat  plate  results  reported  by  Vanier  and  Tien  (1968),  which  show  a  sharp  minimum 
between  5°  <Too<  6°C.  To  ascertain  the  effect  of  sphere  diameter,  Vanier  and  Tien  (1970)  proposed 
a  correlation  of  the  sphere  results  with  those  from  theoretical  analysis  of  a  melting  plate  by 

Nuv/Nu  =C{L/D)3'4  (113) 

where  L  and  D  are  the  characteristic  height  of  the  plate  and  the  diameter  of  the  sphere.  The  least- 
squares-fitted  constant  C  was  found  to  be  1 . 106  ±  0. 1 44.  The  scaled-up  experimental  data  are  shown 
in  Figure  75,  which  clearly  shows  that  the  melting  sphere  behaves  very  similarly  to  a  melting  flat  plate 
and  that  if  all  the  transfer  parameters  are  equal  (including  temperature,  characteristic  length,  and 
surface  area),  about  1 1%  more  heat  is  transferred  to  the  plate  than  to  the  sphere.  This  is  the  effect  of 
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Figure  76.  Variation  of  mean  Nusselt  numbers  with 
ambient  temperature  (experimental  results).  The 
solid  curve  is  from  the  analytical  work  of  Gebhart  and 
Mollendorf  (1978),  the  dashed  curve  is  the  prediction  with 
the  Boussinesq  approximation  (after  Bendell  and  Gebhart 
1976). 

curvature  on  the  flow  velocities  and  is  in  good  agreement  with  the  analytical  results  reported  by  Merk 
and  Prins  (1954)  for  non-melting  free  convection  systems  without  maximum  density  effects,  i.e., 

m7p/m7  =  1.14  (l/d) 3/4 .  di4) 

The  minimum  Nusselt  number  for  spheres  occurs  at  T ^  =  5.35  ±  0.2°C,  as  compared  to  the  value  of 
5.31°C  based  on  Merk’s  (1954)  theoretical  results. 

The  most  recent  study  of  heat  transfer  and  ice  melting  in  ambient  water  near  its  density  extremum 
was  reported  by  Bendell  and  Gebhart  ( 1 976).  In  their  experiment,  a  vertical  ice  slab  (30.3  cm  high,  1 4.8 
cm  wide  and  3  cm  thick  initially)  was  immersed  in  water  at  a  uniform  bulk  ambient  temperature,  T^. 
Figure  76  shows  the  experimental  results  along  with  the  analytical  results  of  Gebhart  and  Mollendorf 
(1978)  and  those  predicted  with  the  Boussinesq  approximation.  Gebhart  and  Mollendorf  s  work  is 
similar  to  that  of  Vanier  and  Tien  (1968)  except  that  it  gives  a  more  accurate  representation  of  the 
density-temperature  relationship  of  water.  As  Vanier  and  Tien  pointed  out,  the  validity  of  the  simplest 
boundary  layer  theory  becomes  questionable  in  the  inversion  region.  However,  beyond  that  region,  the 
experimental  Nusselt  number  values  are  nearly  equal  to  those  predicted  by  theoretical  analysis. 
Bendell  and  Gebhart  ( 1 976)  reported  that  for  a  melting  vertical  ice  surface,  upflow  occurred  when 
<,  5.6°C.  For  roo>5.5°C,  downflow  was  observed  and  was  found  to  be  in  good  agreement  with  earlier 
results.  Bendell  and  Gebhart  (1976)  found  the  minimum  Nusselt  number  for  the  experimental 
temperature  range  2.2°  27^2  25.2°C  occurred  at  7^  =  5.6°C.  In  the  immediate  neighborhood  of  the 
flow  direction  inversion  TIW,  very  slow  flow  exists  with  the  effective  Grashof  number  becoming  very 
small.  Therefore,  the  validity  of  the  simplest  boundary-layer  theory  becomes  questionable,  and  no 
theoretical  results  forthis  regime  were  found  in  the  literature.  After  multiplication  by  a  factor  of  (0. 1 02/ 
0.303 )1/4(sol id  curve  in  Fig.  76),  the  theoretical  results  of  Gebhart  and  Mollendorf  (1978)  generally 
compared  remarkably  well  with  those  reported  by  Vanier  and  Tien  (see  solid  curve  in  Fig.  75),  even 
though  a  rather  elaborate  and  more  accurate  density-temperature  relationship  of  water  was  claimed 
to  be  used  in  Gebhart  and  Mollendorf  s  study. 
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Figure  77.  Velocity  profiles  for  =  2.0  °C  (Wilson  and  Vyas  1979). 


Wilson  and  Vyas  (1979)  reported  the  only  experimental  results  concerning  a  velocity  profile  near 
a  vertical  ice  surface  melting  into  fresh  water  at  temperatures  from  2°  to  7°C.  The  results  suggested 
that  upward  flows  exist  for  water  temperatures  below  4.7°C.  Completely  downward  flowing 
boundary  layers  are  suggested  for  temperatures  about  7°C.  At  4.7°C  an  oscillatory  dual  flow  situ¬ 
ation  begins  to  occur.  As  the  temperature  is  increased,  this  phenomenon  is  increasingly  prevalent, 
with  increasing  downward  velocities  reaching  a  maximum  at  5.6°C.  Figure  77  shows  a  comparison 
between  the  Wilson  and  Vyas  (1979)  experimental  data  and  predictions  obtained  from  free 
convective  boundary  layer  theories  by  Bagley  et  al.  (1972)  and  Vanier  and  Tien  (1970).  The  work 
of  Bagley  et  al.  predicts  a  maximum  velocity  higher  than  that  of  the  data  and  a  significantly  thinner 
boundary  layer  thickness.  Disagreement  may  be  attributed  to  the  high  degree  of  sensitivity  of  the 
analysis  to  the  fluid  density. 

As  indicated  in  Figure  78a  as  the  water  temperature  rises  to  3.2°C,  the  net  buoyancy  forces  in¬ 
crease  resulting  in  an  increase  in  the  maximum  velocity.  At  4.4°C,  the  flow  is  still  upward,  although 
the  maximum  velocity  is  lower  and  the  boundary  layer  is  slightly  thinner  than  at  3.2°C  (Fig.  78b). 
because  of  the  downwards  buoyancy  occurring  in  the  portion  of  the  profile  with  the  temperature 
between  4°  and  4.4°C.  At  4.7°C  (Fig.  78c)  a  similar  but  stronger  tendency  is  noted  as  the  result  of  the 
increased  downward  buoyancy  existing  in  the  outer  portion  of  the  boundary  layer.  The  distinct 
difference  in  the  two  maximum  velocities  suggested  the  oscillating  nature  of  the  water  flow  at  this 
temperature. 

At  5°C  (Fig.  79a)  the  oscillatory  dual  flow  characteristics  are  clearly  established,  and  it  can  be  noted 
that  the  upward  flowing  portion  is  substantially  thinner  and  is  accompanied  by  a  decrease  in  the 
maximum  upward  velocity  when  downflow  exists.  This  trend  is  continued  for  T- 5.3°,  5.6°  and  5.9°C 
as  shown  in  Figures  79b,  c  and  d,  respectively.  In  each  case,  oscillatory  dual  flow  exists  continuously; 
i.e.,  upward  flows  were  not  observed  at  any  time  beyond  the  near-wall  upward  flow  region.  At  5.3°C, 
the  position  of  maximum  downward  velocity  oscillates  between  7.3  and  1 1 .5  mm  from  the  wall.  For 
5.6°C,  it  lies  between  6  mm  and  18  mm,  but  at  5.9°C,  this  variation  is  limited  to  5.2  to  8.0  mm.  In 
addition,  the  profiles  at  5.9°C  (Fig.  79d)  pose  much  greater  self-similarity  than  those  of  5.3°  and 
5.6°C.  Thus  the  flow  phenomena  for  5.3°  and  5.6°C  seem  to  be  substantially  more  complex  than  that 
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Figure  78.  Velocity  profiles for  a)  T..  =  3.2  V,  b)  =  4.4  V  andT^ 

4.7V  (Wilson  and  Vyas  1979). 


at  5.9°C.  In  both  Figures  79b  and  79d  there  are  small  regions  of  the  order  of  1  to  2  mm  in  thickness 
near  the  wall  where  upwards  flow  exists,  while  at  5.6°C  much  smaller  velocities  were  observed 
within  the  same  region  but  with  much  lower  maximum  velocity.  Figure  79e  shows  the  characteristics 
of  a  relatively  stable  downward  flow  regime.  In  this  case,  the  portion  of  the  boundary  layer  with 
upward  buoyancy  is  approximately  0.5  to  1 .0  mm  in  thickness.  The  downward  buoyancy  forces  in 
the  outer  portion  are  much  stronger  and  the  whole  flow  is  downward,  because  the  viscous  forces, 
which  tend  to  draw  the  fluid  downward,  approximately  balance  the  upward  buoyancy  forces  near 
the  ice.  The  velocity  gradient  is  approximately  zero  at  the  ice  surface. 
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DISCUSSION  AND  CONCLUSIONS 

This  review  covers  only  the  problems  associated  with  the  anomalous  density-temperature 
relationship  of  water  contained  in  various  geometrical  systems.  The  discussion  and  conclusion  are 
classified  into  two  subjects:  1)  onset  of  convection,  2)  flow  patterns  and  natural  convective  heat 
transfer. 

Onset  of  convection 

The  criterion  for  the  onset  of  convection  of  a  confined  horizontal  layer  that  contains  a  density 
extremum  was  found  both  experimentally  ( Y en  1 968,  Yen  and  Galea  1 969)  and  analytically  ( Veronis 
1963,  Sun  et  al.  1969,  Merker  et  al.  1979)  to  not  be  a  constant  value,  as  in  the  classical  Benard 
problem,  but  dependent  on  the  thermal  boundary  conditions.*  This  is  evident  in  Figure  4,  in  which 
Rac  was  plotted  explicitly  vs  T]  or  T2  In  the  case  of  melting  from  the  top,  the  higher  the  values  of 
Tv  the  greater  RaQ  becomes;  in  other  words,  the  farther  removed  temperature  T2  is  from  4°C,  the  less 
prone  the  layer  is  to  begin  the  onset  of  convection.  In  the  case  of  melting  from  below,  as  temperature 
T,  increases,  RaQ  decreases  exponentially  and  asymptotically  approaches  the  value  <=  1708,  as  re¬ 
ported  in  the  classical  Benard  problem.  This  is  evident  because,  as  f,  becomes  higher  and  higher, 
the  buoyancy  forces  created  by  temperature  difference  AT  (=  T,  -  7"max)  possess  a  much  stronger 
influence  on  the  layer  stability  than  the  effect  produced  by  the  density  extremum  (i.e.,  at *»  4°C),  and 
subsequently  the  continuously  forming  layer  behaves  like  a  normal  fluid  (i.e.,  there  is  a  monotonic 


*  If  the  layer  is  formed  by  phase  transition,  then  one  of  the  thermal  boundaries  is  at  the  ice  melting  point,  i.e., 

T  =0°C. 
m 
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density-temperature  relationship)  as  in  the  Benaid  problem.  On  the  other  hand,  as  Tx  decreases  and 
approaches  the  temperature  of  the  density  extremum,  Rac  increases  and  approaches  the  limiting 
value  of  infinity .  This  is  also  expected,  since  if  Tj  is  maintained  at  *  4°C,  the  water  has  its  greatest 
density  at  the  lower  boundary  and  the  water  layer  will  always  remain  stable. 

In  the  case  of  a  water  layer  formed  by  melting  ice  from  above,  the  trend  of  variation  of  Rac  with 
boundary  temperature  is  simply  reversed.  The  higher  the  temperature  T2,  the  greater  Rac  becomes. 
This  can  be  explained  by  noting  that  if  T2  is  maintained  in  the  range  of  0  <  r224°C,  the  entire  layer 
is  unstable  because  the  higher  density  water  will  overlie  the  less  dense  water  underneath,  and  will 
consequently  result  in  lower  RaQ  values.  If  T2  is  maintained  at  a  higher  temperature  than  4°C,  only 
a  fraction  of  the  layer  (=  4 H/T2)  is  potentially  unstable,  and  the  layer  is  less  prone  to  onset  of 
convection.  It  seems  that  the  RaQ  value  grows  increasingly  greater  as  the  effect  of  the  density 
extremum  becomes  less  pronounced.  It  is  also  interesting  to  note  that  the  two  Rac  curves  intersect 
at  exactly  T{  =  T2 = 8°C,  which  clearly  indicates  that  under  these  particular  thermal  conditions,  the 
two  systems  are  identical  and  have  a  unique  RaQ  value  regardless  of  how  the  water  layer  was  formed. 

In  the  case  of  a  confined  pure  water  layer,  Merker  et  al .( 1 979)  reported  Rac  values  for  both  constant 
temperature  and  constant  heat  flux  thermal  boundaries.  The  effect  of  Tt  and  T2  on  Rac  seems  to  be 
similar  but  with  higher  values  of  Rac  for  the  cases  of  Tw  =  constant  (see  Fig.  6  and  7).  They  indicated 
that  the  values  for  Rac  depend  on  the  specific  values  off,  and  T2.  A  few  comparisons  between  the 
results  from  Sun  et  al.  and  Merker  et  al.  were  found  to  be  in  good  agreement. 

The  stability  problem,  equivalent  to  the  case  of  melting  from  above,  was  studied  by  Seki  et  al. 
( 1 977)  with  the  added  condition  of  a  free  surface  maintained  either  at  T2  >  4°C  or  T2  <,  4°C  and  T{ 
at  0°C.  For  T2  <  8°C,  they  demonstrated  both  analytically  and  experimentally  that  the  criterion  of 
hydrodynamic  stability  in  a  water  with  a  density  inversion  is  dependent  on  the  free  water-surface 
temperature  T2and  increases  as  surface  tension  increases  and  decreases  as  T2  is  lowered,  resulting 
in  a  reduction  of  the  stable  layer  thicknesses.  On  the  other  hand,  for  T2  S  8°C,  the  criterion  is  found 
to  be  independent  of  T-,  and  approaches  asymptotically  a  limiting  value  of  a  modified  critical 
Rayleigh  number,  i.e.,  Rac  =■  500,  as  predicted  by  Sun  et  al.  (1969). 

Hassab  and  Sorour  ( 1982)  reported  the  first  analytical  study  on  the  stability  of  the  conduction 
regime  of  natural  convection  in  a  vertical  melt  layer  formed  by  melting  ice.  Their  stability  criterion 
is  expressed  in  terms  of  critical  Grashof  number  Grc.  The  values  of  Gr,  (and  the  corresponding  critical 
melting  thickness  hc )  are  dependent  upon  the  stepped  wall  temperature  Tx,  such  that  as  T,  is  in¬ 
creased  the  change  in  heat  transfer  mode  from  conduction  to  convection  is  enhanced  for  ?11  values 
of  T,  in  the  range  of  1  °  to  30°C:  They  found  that  the  instability  sets  in  as  vertical  traveling  waves, 
with  the  secondary  flow  occurring  as  two-column  waves  for  T,  <  7. 1  °C  and  Tx  >  9.4°C,  and  as  three- 
column  waves  for  7. 1°C  <  T,  <  9.4°C. 

Flow  patterns  and 

natural  convective  heat  transfer 

The  most  striking  phenomenon  that  results  from  the  presence  of  a  fluid  density  maximum  within 
a  confined  area  is  the  creation  of  an  unusual  temperature  distribution  and  its  associated  cellular 
formation  and  flow,  which  are  directly  related  to  geometrical  arrangement  (i.e.,  in  a  confined 
horizontal  layer,  in  a  vertical  gap,  in  a  cylindrical  annulus,  or  in  a  rectangular  and  square  enclosure). 
In  a  horizontal  layer  formed  by  melting  ice  (in  this  case,  one  in  which  the  boundary  is  ice  at  its 
melting  temperature),  the  transient  temperature  and  its  pseudo-steady  temperature  distribution  can 
be  represented  as  shown  in  Figures  67  and  68.  It  can  be  noted  that  the  unique  features  of  these 
distributions  are  the  formation  and  expansion  of  the  constant  temperature  region.  The  heat  flux 
across  the  water/ice  interface  was  found  to  be  a  weak  function  of  temperature  for  the  melt  layer  for 
the  case  of  melting  from  above  and  a  linear  function  of  temperature  for  the  case  of  melting  from 
below. 

Lankford  and  Bejan  ( 1 986)  developed  a  unique  heat  transfer  correlation  based  on  scale  analysis 
for  water  near  4°C.  The  curve  shown  in  Figure  27  was  obtained  from  eq  35,with  C,  =  0.31  and  C2 
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=  O.S,  and  a  set  of  specific  Rayleigh  numbers  based  on  the  hot  and  cold  sides  boundary  layers, 
respectively.  Without  this  modification  the  ordinary  correlation  for  normal  fluids  (or  water  warmer 
than  4°C),  e.g.,eq25  or  Figure  25,  failed  to  represent  the  data,  including  the  results  for  temperatures 
near  4°C  or  below. 

Figures  29a,  b,  c,  and  d  typically  demonstrated  the  temperature  distribution  and  cell  rotation 
pattern  in  a  square  enclosure  (aspect  ratio,  A  =  1,  K  =  rjr  =  1;  no  curvature  effect)  as  function  of 
values  of  R\  when  R' = 0.5  (i.e.,  when  the  maximum  density  is  in  the  middle  plane).  The  temperature 
distribution  as  well  as  the  cell  size  are  nearly  identical  except  rotating  in  the  opposite  direction. 
When  R7  =  1.0,  (i.e.,  the  hot  wall  temperature  is  equal  to  the  Tm),  the  water  behaves  as  a  normal 
Boussinesq  fluid  except  for  the  direction  of  circulation.  The  effect  of  A  on  temperature  distribution 
and  flow  pattern  is  clearly  demonstrated  in  Figure  31 .  As  A  increases,  the  convection  heat  transfer 
occurs  only  at  the  top  of  the  annulus,  and  not  at  the  near-bottom,  because  of  the  dual  cell  structure. 
The  effect  of  K  and  Ra'  on  the  heat  transfer  at  the  inner  wall  (relative  to  conduction)  is  shown  in 
Figure  33.  From  this  study,  it  is  evident  that  density  inversion  phenomena  are  altered  substantially 
by  the  curvature  of  the  annulus.  A  particular  steady  flow  structure  is  determined  by  the  combination 
of  the  inversion  parameter  R'  and  the  curvature  K.  A  transition  from  inner  to  outer  convective  cell 
dominance  can  be  accomplished  by  either  increasing  R'  or  by  increasing  the  curvature  K  of  the  annulus. 

Seki  et  al.  ( 1 975)  described  the  only  experimental  study  of  the  effect  of  water  near  4°C  on  natural 
convective  heat  transfer  between  two  horizontal  concentric  cylinders.  They  reported  that  the  effect 
of  density  inversion  is  unexpectedly  large  and,  especially  when  two  counter-eddies  of  approxi¬ 
mately  the  same  size  coexist  in  the  gap,  the  average  Nusselt  number  indicated  a  minimum  value  and 
the  minimum  Nu  increasing  as  the  gap  width  increased  (see  Fig.  39).  However,  there  seemed  to  be 
no  definite  pattern  as  to  the  temperatures  where  the  minimum  heat  transfer  occurred. 

The  analytical  work  of  Nguyen  et  al.  (1982)  confirmed  the  experimental  findings  of  Seki  et  al. 
( 1 975).  The  pattern  of  circulation  and  angular  velocity  are  functions  of  inversion  parameters  y '  as 
well  as  outer  to  inner  radius  ratio,  i.e.,  R '  and  the  Rayleigh  number.  A  slight  variation  in  y' creates 
a  radical  change  in  cellular  configuration.  The  minimum  coefficient  of  convective  heat  transfer  is 
found  to  deviate  from  y '  =  -1  as  R'  increases  (see  Figs.  40, 41 , 42  and  43). 

The  extension  of  this  woik  by  Vasseuret  al.  (1983)  further  delineated  the  temperature  distribution 
and  the  circulation  pattern  within  the  cylindrical  annulus ,  especially  for  high  Rayleigh  numbers.  For 
R'  =  2.6  and  Ran  varying  from  2  x  103  to  7  x  104,  an  overall  minimum  Nu  was  found  to  occur  at  y ' 
=  -0.85.  On  the  other  hand  for  Ran  =  1  x  104,  the  minimum  Nu  occurs  at  y  *  >  -1  but  approaches 
y'  =  -1  as  R'  increases  from  1 .75  to  2.6. 

The  work  of  Watson  ( 1 972)  and  later  Lin  and  Nansteel  ( 1 987a)  on  the  effect  of  density  inversion 
on  the  temperature  distribution  and  heat  transfer  in  a  square  enclosure  is  a  special  case  of  that  for 
rectangular  enclosures.  Watson  demonstrated  that  the  variation  of  Nu  with  the  warm  temperature 
boundary  7”,  between  the  full  equation  and  constant  viscosity  model  are  negligible  for  T,  <  8°C  but 
deviate  from  each  other  as  T,  >  ~  9°C.  On  the  other  hand,  the  Boussinesq  model  displayed  a 
completely  different  heat  transfer  phenomenon  (see  Fig.  47).  For  a  square  enclosure,  for/?'  values 
ranging  from  0  to  1,  the  average  Nu  was  found  to  be  symmetrical  with  respect  to  R'  =  0.5,  with  a 
minimum  at  R'  =  0.5,  and  increases  as  Ra  increases. 

For  rectangular  enclosures,  Desai  and  Forbes  (1981)  numerically  calculated  that  the  Nusselt 
number  is  always  higher  for  the  T-T^  range  of  0°  to  8°C  than  for  the  T-T^  range  of  2°  to  6°C 
(although  for  both  cases  R'  =  1/2)  independent  of  the  temperature-density  representation  and  aspect 
ratio.  Nansteel  et  al.  ( 1 987)  also  found  that  the  value  of Nu  is  very  sensitive  toR'  as  well  as  the  aspect 
ratio.  The  minimum  Nu  occurs  at  R'  =  1/2,  and  it  increases  with  increasing  A. 

The  effects  of  4°C  on  the  natural  convective  heat  transfer  and  temperature  distribution  with  initial 
temperatures  at  4°  and  8°C  were  reported  by  Forbes  and  Cooper  (1975)  who  cooled  water  from  the 
top  with  either  a  rigid  boundary  condition  at  constant  temperature  or  a  free  water-air  surface  with 
constant  convective  heat  transfer  coefficient.  They  demonstrated  that  the  eddy’ s  formation  is  closely 
related  to  the  position  of  the  maximum  isotherm.  At  the  beginning,  one  cell  was  formed  below  the 
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4°C  line  (square  enclosure).  As  the  4°C  curve  descended  farther,  two  cells  were  formed  on  top  of 
each  other,  with  part  of  the  4°C  line  as  the  boundary  and,  finally,  the  two  cells  returned  to  one  eddy 
(or  one  cell)  as  the  4°C  line  progressed  through  the  bottom  of  the  container.  For  W/D  >  1 ,  the  same 
mechanism  of  cell  formation  was  observed,  but  the  number  of  eddies  formed  was  not  in  any  way 
related  to  the  values  of  W/D.  (Based  on  limited  data,  for  W/D  =  3,  three  eddies  with  unequal  size  were 
formed;  however,  for  W/D  =  6,  seven  eddies  were  noted.) 

Yen  (1968)  and  Yen  and  Galea  (1969)  performed  the  only  experimental  work  on  temperature 
distribution  and  heat  transfer  in  a  melt  layer  formed  by  melting  ice  (with  ice  initially  at  its  fusion 
temperature).  Four  cases  of  melting  layers  were  set  up,  with  T,  S4°C  (melting  from  below)  and  T2 
S4°C  (melting  from  above).  In  the  case  off,  =  4°C  the  whole  melt  layer  was  always  stable,  but  for 
T2  =  4°C  the  whole  layer  was  unstable.  For  7^  and  T2  >  4°C,  the  melt  layer  consisted  of  one  stable 
and  one  unstable  region.  The  effects  of  these  unstable  regions  on  temperature  distribution,  heat 
transfer  and  onset  of  convection  were  found  dependent  on  mode  of  melting.  The  most  striking 
feature  in  the  ice-melting  system  was  the  formation  andexpansion  of  the  constant  temperature  zone, 
which  had  a  temperature  of  about  3.2°C,  regardless  of  the  value  o  fTr  but  had  a  dependency  on 
in  the  case  of  melting  from  below. 

For  the  case  of  ice  melting  in  bulk  water,  experimental  and  analytical  studies  reveal  the  existence 
of  split-flow  at  the  inversion  temperature  Tiv,  (i.e.,  for  7^  <  Tjv  the  boundary  layer  flow  is  upward, 
and  for  >  Tjv  the  entire  boundary  layer  is  downward).  The  inversion  temperature  was  found  to 
be  about  5°  or  6°C  for  a  great  number  of  phase-transition  geometries,  but  not  at  the  temperature  of 
maximum  density  (=  4°C). 
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